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TRANSIENT  FIELDS  PRODUCED  IN  HETEROGENEOUS  BODIES 
BY  ELECTROMAGNETIC  PULSES 

STATEMENT  OF  THE  PROBLEM 

This  effort  examines  the  question  of  electromagnetic  field  structure 
induced  in  heterogeneous  bodies  by  pulse  electromagnetic  radiation.  Of  spe¬ 
cial  interest  is  the  case  of  a  sinusoidally  varying  incident  field  of  finite 
duration. 

To  determine  the  structure  of  induced  electromagnetic  fields,  the  problem 
was  modeled  mathematical ly.  The  assumptions  made  and  the  resulting  model  are 
given  below.  As  the  model  is  described,  the  nature  of  the  media  and  incident 
fields  to  be  considered  will  be  made  clear. 

The  medium  considered  in  this  work  is  one  dimensional  and  of  finite  depth 
L.  Thus,  it  is  assumed  to  vary  spatially  only  in  the  z  direction,  extending 
from  z= 0  to  z= L.  The  permittivity,  e(z),  and  conductivity,  o(z),  of  the 
medium  are  piecewise  continuous  functions  as  pictured,  for  example,  in  Fig¬ 
ure  1.  By  allowing  e  and  a  to  be  piecewise  continuous,  a  layered  medium  can 
be  considered.  Free  space  is  assumed  to  exist  outside  of  the  medium,  so 
e ( z ) =£ g  and  o(z)=0  for  z<0  and  z>L.  Finally,  the  medium  is  assumed  to  be  non- 
dispersive  in  that  the  permittivity  and  conductivity  are  independent  of  the 
induced  field. 

An  electromagnetic  wave  propagating  along  the  z-axis  normal  to  the  medium 
has  transverse  electric  field,  E(z,t),  satisfying 

Ezz"e^MoEtt‘a^MoEt=0,  -»<z<~,  -»<t<«  (1) 

which  follows  from  Maxwell's  equations.  Here,  t  denotes  time,  is  the  con¬ 
stant  magnetic  permeability,  and  subscripts  denote  partial  derivatives.  The 
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Figure  1.  Example  of  permittivity  and  conductivity  profiles 
for  a  one-dimensional,  three-layer  medium. 


incident  electromagnetic  pulse,  E1,  is  assumed  to  propagate  from  left  to  right 
and  first  impinges  on  the  medium  at  time  t=0.  Thus,  the  initial  condition  for 
equation  (1)  is 

E(z,t)=Ei(z-ct)  for  t<0  (2) 

where  E*  (z-ct)=0  for  z-ct>0  and  c  is  the  speed  of  light  in  free  space.  Fig¬ 
ures  2  and  3  show  what  some  typical  incident  pulses  may  look  like. 


The  problem  considered,  then,  is  that  of  determini nq  the  solution  E(z,t) 
of  the  system  of  equations  (1)  and  (2)  for  0<z<L,  t>0  for  any  given  E’fz-ct), 
e(z),  o(z). 


SOLUTION  TECHNIQUE 

The  problem  modeled  above  has  been  solved  numerically  by  the  HATS  compu¬ 
ter  program  (High  Accuracy  Temporal  Solution).  A  detailed  description  of  the 
mathematical  foundations  of  the  numerical  methods  used  in  the  program  can  be 
found  in  reference  8,  and  preliminary  results  are  given  in  references  4-7.  A 
brief  summary  of  reference  8  is  given  in  Appendix  A.  A  description  of  the  use 
and  capabilities  of  this  program  is  given  in  Appendix  B,  which  is  a  self- 
contained  user's  manual  for  HATS. 

Note  that  the  numerical  solution  is  based  on  exact  analytical  represen¬ 
tations  of  the  fields  involved,  with  no  recourse  to  low  frequency  approxima¬ 
tions.  The  technique  used  has  been  found  to  be  highly  stable. 
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SAMPLE  PROBLEM 

As  an  example  of  the  results  that  can  be  obtained  via  the  HATS  program, 
consider  the  three-layer  medium  shown  in  Figure  4,  where  it  is  first  assumed 
that  o(z)=0.  A  1000-MHz  incident  field  of  20-ns  duration  impinges  on  the 
medium  from  the  left  at  time  t=0.  Thus,  at  some  time  t<0  the  incident  field 
looks  like  that  shown  in  Figure  5. 

The  resulting  field  inside  the  medium  was  generated  by  the  HATS  program 
for  the  time  interval  t=0  to  t=25  ns.  Figures  6,  7,  and  8  show  the  internal 
fields  as  a  function  of  time  at  three  specific  locations  (z=.714  cm,  z=1.067 
cm,  z=2.100  cm)  in  the  medium.  The  locations  z=.714  cm  and  z=2.100  cm  (Fiqs. 
6  and  8)  were  chosen  because  they  represent  in  some  sense  the  "worst  case" 
behavior,  in  that  the  transients  at  these  points  had  the  qreatest  amplitude 
relative  to  steady  state  amplitude  of  any  points  in  the  medium.  Flqure  7 
(z=1.067  cm)  represents  the  "best  case"  behavior  in  that  the  ratio  of  tran¬ 
sient  amplitude  to  steady  amplitude  was  the  smallest.  At  all  points  in  the 
medium,  transient  amplitude  was  greater  than  steady  state  amplitude. 
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Figure  4.  Permittivity  profile  for  the  sample  problem 


20  cycles 


f  fc<o 


Figure  6.  Internal  field  (z=0.714  cm)  vs.  time  (ns),  zero  conducti vity. 


time  (ns),  zero  conduct 
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Now  consider  a  medium  with  permittivity  profile  given  in  Figure  5  and 
conductivity  profile  given  in  Figure  9.  The  fields  at  z=.714  cm,  z=1.067  cm, 
z=2.100  cm  (Figs.  10-12)  fairly  accurately  give  worst  and  best  behavior. 
Again,  transient  amplitude  was  everywhere  greater  than  steady  amplitude, 
although  the  effect  of  nonzero  conductivity  reduces  the  ratio  of  transient  to 
steady  amplitude. 


<rc*) 


Figure  9.  Conductivity  profile  for  the  sample  problem. 


Additional  information  regarding  these  plots  (how  the  specific  values  of 
z  were  chosen,  views  of  the  entire  field  in  the  medium  at  particular  times, 
etc.)  and  a  detailed  discussion  of  how  they  were  obtained  are  given  in  Appen¬ 
dix  B.  The  problem  studied  in  this  section  is  only  a  sample  problem;  any  per¬ 
mittivity  and  conductivity  profiles  and  incident  field  can  be  input  into  the 
HATS  program,  modulo  the  restrictions  called  out  in  Appendix  B. 
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CONCLUSIONS 


Use  of  the  HATS  computer  program  has  shown  that  relatively  large  ampli¬ 
tude  transients  can  occur  in  a  medium  when  a  transient  signal,  such  as  that 
from  a  phased  array  system,  impinges  on  the  medium.  Although  the  model  used 
to  determine  this  is  a  fairly  simple  one,  the  pattern  shown  here  should  extend 
to  more  realistic  situations.  Nonetheless,  further  work  is  needed  to  defi¬ 
nitely  establish  the  presence  and  magnitude  of  these  transients  for  other 
geometries  and  for  the  important  case  of  a  frequency  dependent  medium,  such  as 
human  tissue. 
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APPENDIX  A.  OUTLINE  OF  SOLUTION  TECHNIQUE 

The  problem  considered  In  this  work  is  modeled  by 

Ezz-e(z)u0Ett"a(z)woEta°*  -<*<-.  -°°<t<00*  (A.l) 

It  is  assumed  that  e(z)  and  0(2)  are  given  and  that  e(z)=e0>  o(z)=0  for  z<0 
and  z>i.  Finally ,  the  solution  of  equation  (A.l)  is  given  by 

E(z,t)=E1  (z-ct)  f2  '<0  (A. 2) 

where  E1  (z-ct)=0  for  7,-ct>0.  The  solution  of  the  initial  value  problem 
(A.l,  A. 2)  must  be  found  for  0£z<L  and  t>0. 

To  facilitate  the  numerical  solution  of  equation  (A.l),  the  transforma¬ 
tions 


x=x(z)=(/*{e(s)u{)}ids)/Jl  (A. 3) 

T=t /l  (A. 4) 

u(x,t  )=E(z,t)  (A. 5) 

are  used  to  convert  equation  (A.l)  to 

UXX‘Utt+A(x)ux+B(x)ut=0,  -»<X<»,  (A. 6) 
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where 


A(x)=-^{e(z)M0f* 

B(x)=-ia(z)/e (z) 

and  i  is  such  that  x(L)=l.  In  terms  of  this  transformed  problem,  the  solution 
u(x,t)  of  equation  (A. 6)  must  be  found  for  0<x<l,  t>0  given  that 

u(x,t)=u^ (x-t)  for  t<0  (A. 7) 

where  u1(t)=E1(cHt). 

The  advantage  of  the  change  of  variables  (A.3)-(A.5)  is  that  the  charac¬ 
teristics  of  equation  (A. 6)  are  straight  lines.  On  the  other  hand,  while  the 
solution  E  of  equation  (A.l)  has  continuous  first  derivatives,  the  solution  u 
of  equation  (A. 6)  has  discontinuities  in  its  x  derivative.  To  see  this, 
assume  that  e(z)  is  discontinuous  at  the  points  z=z.j,  i=0,  1»  •••»  NL  where 

0=zq<z1<z2<...<z^l=L.  Let  x.j=x(Zj)  for  i=0,  1 . NL.  Then  the  change  of 

variable  (A. 3)  implies  that 

ciVV’  T )=ux^xi ^  (A-8) 

for  i=0,  1,  ....  NL  where 

Ci-MvlM*,-)!*- 

One  technique  for  solving  the  problem  (A.6)-(A.8)  is  by  using  finite  dif¬ 
ferences  [1].  However,  in  using  an  explicit  differencing  scheme  on  this  prob¬ 
lem,  roundoff  error  accumulated  quickly  and  prevented  accurate  solutions  for 
incident  pulses  of  realistic  duration.  An  alternate  approach  is  to  use  the 
time-domain  integral  equation  technique  of  Bolomey  et  al.  [3]  with  suitable 
modifications  to  take  into  account  the  finite  depth  of  the  medium.  However, 


this  technique  proved  to  be  susceptible  to  numerical  dispersion  in  that  trun¬ 
cation  error  propagates  through  the  solution  at  nonphysical  speeds,  causinq 
spreading  of  the  fields.  Again,  for  pulses  of  realistic  duration,  this  cre¬ 
ates  highly  unreliable  solutions. 

The  above  two  solution  techniques  are  based  on  the  familiar  space-time 
diagram  shown  in  Figure  A-l,  which  displays  the  fact  that  the  solution  u  at 
the  point  (x,t)  depends  on  the  solution  in  a  certain  interval  (the  domain  of 
dependence)  on  the  x-axis  determined  by  the  characteristic  lines  emanatinq 
from  (x,t).  Note,  however,  that  an  equally  valid  solution  technique  is 
afforded  through  the  space-time  diagram  of  Figure  A-2.  To  understand  the  data 
required  for  the  formulation  of  the  problem  suggested  by  Figure  A-2,  note  that 
the  incident  electromagnetic  field  E1  (z-ct)  gives  rise  to  a  transmitted  field 
E^z-ct)  for  z> L.  In  the  transformed  problem  (A.6)-(A.8),  the  incident  plane 
wave  u1  (t)=e’ (cat)  produces  a  transmitted  plane  wave  ut(x-i)  where  u^x) 
=E*  (L-ct+c2x)  and  u1 (x)=ut  (t)=0  for  x>0.  Since  u(l,x)-ut  (1-t),  it  follows 
that  the  data  required  for  the  formulation  of  the  problem  shown  in  Figure  A-2 
is  u1  (x)  for  x<0.  The  advantage  of  this  formulation  over  that  suggested  by 
Figure  A-l  is  that  exact  analytical  expressions  exist  for  the  solution  u(x,x) 
in  terms  of  the  data  u^  and  these  expressions  give  rise  to  highly  accurate 
numerical  schemes. 

Pursuing  the  approach  suggested  by  Figure  A-2,  it  is  first  necessary  to 
determine  the  transmitted  field  u^  corresponding  to  a  given  incident  field  u  . 
In  reference  5  these  fields  are  shown  to  be  related  by  the  equation 


The  constants  Wj  and  rj  are  known  and  are  determined  by  the  properties  of  the 
medium  (formulas  are  given  in  reference  5).  The  kernel  function,  W,  also 
depends  only  on  the  properties  of  the  medium  and  is  Independent  of  the  fields 
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Assuming  for  the  moment  that  W  is  known,  equation  (A. 9)  is  a  delay 
Volterra  integral  equation  for  u*.  Thus,  it  is  equivalent  to  a  Volterra  inte¬ 
gral  equation  of  the  second  kind  (as  shown  in  reference  6),  so  stable  numeri¬ 
cal  techniques  for  solvlnq  equation  (A. 9)  for  u*  are  available  (see  refer¬ 
ence  2).  Consequently,  the  data  required  for  the  formulation  in  Figure  A-2 
can  be  considered  known. 

Reference  8  shows  that  for  any  x  in  the  interval  [ xNt__ j  ,  l]  and  for  any 
t,  the  solution  u(x,t)  of  equation  (A. 6)  is  given  by 

2u(x,T)=exp(l/2/^(A(s)-B(s))ds}{(cNL+l)ut(x-T) 

-(cNL-l)ut(2-x-T)exp(/^B(s)ds} 


+J2’V(y-T)N(x,y,l,)dy}.  (A. 10) 

x 

Here,  N(x,.y,l)  is  the  solution  of  the  characteristic  initial  value  problem 

Nxx-Nyy+B(x)(Nx+Ny)+D+^x>N=0*  xNL-lixixNL=1’  (A.ll) 

2N(x,x,l)=cNLB(l)-A(l)-(cNL+l)A+(s)ds,  (A. 12) 

2N(x,2-x,l)=(cNLB(l)-A(l)+(cNL-l)/V(s)ds)  exp{/*B(s)ds}  (A. 13) 


where 

A ( 1 ) =A ( 1 - ) ,  B( 1 )=B( 1- ) 
and 

0±=l/4(B2-A2)+l/2(-A'±B'). 
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The  problem  (A. 11 )-(A.13)  can  be  easily  solved  numerically  by  the  method 
of  characteristics  [1].  Thus,  u(x,t)  can  be  determined  from  equation  (A. 10) 
for  any  x  in  [xNL_j,  l]  by  means  of  a  simple  quadrature.  To  determine 
u(x,t)  for  x  in  [x^.^  »  XNL-1  1*  re^erence  8  shows  that  a  system  of  formulas 
analogous  to  equations  (A.10)-(A.13)  utilizes  the  data  u( xnl_i »  T)*  Continuing 
in  this  manner,  quadrature  formulas  can  be  used  to  determine  u(x,t)  for  any  x 
in  the  interval  [0,1].  In  the  program  HATS,  the  trapezoidal  rule  is  used  to 
perform  these  quadratures. 

The  remaining  point  to  be  considered  is  determination  of  the  kernel  func¬ 
tion,  W,  which  appears  in  equation  (A. 9).  If  the  transmitted  field  is  a  Dirac 
delta  function, 

ut(T)=6(x)  (A. 14) 

then,  according  to  equation  (A. 9),  the  "incident  field"  that  would  qenerate 
such  a  transmitted  field  is  the  kernel  W,  modulo  delta  function  singulari¬ 
ties.  But  if  the  transmitted  field  is  assumed  to  be  given  by  equation  (A. 14), 
then  the  technique  outlined  above  enables  one  to  determine  the  corresponding 
field  u(0,t)  for  all  r.  From  this  it  is  possible  to  separate  out  the  incident 
field  that  generates  (A. 14)  and  hence  obtain  W.  Again,  details  are  given  in 
reference  8. 
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I.  INTRODUCTION 


The  following  general  description  is  meant  to  be  a  user's  guide  to  the 
program  HATS  (High  Accuracy  Temporal  Solution)  and  does  not  dwell  on  mathe¬ 
matical  details  other  than  what  is  needed  to  properly  use  the  program. 

Briefly  stated,  this  program  allows  a  user  to  cause  an  incident  plane 
electromagnetic  wave  to  impinge  normally  on  a  medium  of  finite  depth  and  to 
determine  the  resulting  scattered  and  internal  fields.  The  medium  is  assumed 
to  be  one-dimensional  (i.e.,  stratified  in  the  z  direction)  and  "linear"  in 
that  its  permittivity,  e(z),  and  conductivity,  a(z),  are  independent  of  fre¬ 
quency.  Free  space  is  assumed  to  exist  on  either  side  of  the  medium.  The 
functions  e(z)  and  o(z)  are  to  be  piecewise  continuously  differentiable  and, 
hence,  may  have  jump  discontinuities.  If  a  layer  is  defined  to  be  an  interval 
in  the  medium  on  which  e  and  o  are  continuously  differentiable,  then  the  pro¬ 
gram  HATS  is  capable  of  handling  a  five-layer  medium.  A  two-layer  medium  is 
pictured  in  Figure  B-l. 

The  incident  field  can  be  of  any  form  the  user  desires  with  the  following 
two  provisions: 

1.  The  incident  field  must  be  a  continuous  function. 

2.  The  incident  field  does  not  contact  the  mediun  with  time  t=0. 

Thus,  transient  as  well  as  steady  state  fields  can  be  determined  by 
this  program. 

See  Figure  B-2  for  a  picture  of  some  typical  Incident  fields. 

Details  regarding  proper  Input  Into  the  program  and  form  of  the  output 
will  be  given  in  the  following  sections. 
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Figure  B-2. 


Some  typical  Incident  fields,  E^(i.  -  t),  where  c  Is  the  speed  of 

c 

light  In  free  space. 
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II.  INPUTTING  e  AND  o 


The  program  presumes  that  a  Liouville  transform  has  been  made  on  the 
spatial  variable  z  and  that  a  normalization  has  been  performed  on  the  time 
variable  t.  Therefore,  a  certain  amount  of  massaqinq  must  be  done  before  data 
and  functions  are  put  into  the  program,  and  some  interpretation  of  the  program 
results  is  necessary.  All  data  are  initially  expressed  in  rationalized  MKS 
units. 

Consider  a  medium  of  NL  layers,  where  KNLX5.  Layer  i  extends  from 
z=z.  i  to  1-1. . ,  where  zQ=0  and  ZNL=L.  Assume  that  the  permittivity  and  con¬ 
ductivity  are  given  on  layer  i  by  the  functions  (z)  and  o^(z),  z^_j<z<zt- , 
and  let 


e(z)=e^(z)  for  z^_j<z<z^,  1*1,  •••,  NL. 
The  Liouville  transform 


x=x(z)=(J^{e(s)u0}*ds)/* 


(B.l) 


where 


t=/^|e(s)u0}ids 
and  the  normalization 

T  =  t/£ 

is  used  to  produce  a  problem  in  (x,t)  coordinates  in  which  the  velocity  of 
propagation  is  unity. 

In  the  (x,t)  coordinate  system  the  scattering  medium  extends  from  x=0  to 
x=l  and  is  again  divided  Into  NL  layers.  Properties  of  each  of  these  layers 
(in  terms  of  the  x  variable)  are  put  into  the  program  as  function  statements 
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appearing  In  subroutine  EVAL.  For  each  layer,  four  functions  must  be  speci¬ 
fied.  For  layer  I  (KKNL)  these  functions  are 


AI(x)— ii-{e  (z)u0}  *,  z  <z<z 
dz  I  u  1-1  I 

AIP(x)=  fL  A I ( x ) 
dx 

BI(x)=-toI(z)/eI(z),  Zj_j<z<Zj 

B IP ( x )=  ^_BI(x). 
dx 

These  functions  are  to  be  entered  In  the  appropriate  lines  In  subroutine 
EVAL.  If  NL  Is  less  than  five,  the  statement  functions  which  appear  In  EVAL 
for  layers  NL+1,  •••,  5  are  not  used.  These  unused  statements  may  be  left  In 
the  program. 

As  an  example  of  how  to  determine  these  Input  functions,  let  NL=2,  zQ=0, 

ZjS.Ol  m,  z2=.03  m,  and 


o1(z)=s1  mho/m 
o2(z)=s2  mho/m 

where  e1#  e2,  s2  are  constants  and  cQ  Is  the  permittivity  of  free  space, 
e0*8.854.1o'12  coul2/nm2. 

Then 


KVo*  ds  +  / 


•03,  ,J 

ds 


(O.Ol/ej"  +  0.02/ip/c 
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where  c  Is  the  speed  of  light  In  free  space.  It  follows  that  for  0<z<0.01  m, 

x=x(z)=(z/ep/(tc) 
and  for  0.01<z<0.03, 

x-x(z)=(0.0l/ej*+(z-0.01)/e^)/(tc) . 


Thus, 


Al(x)=0 

AlP(x)=0 

Bl(x)— 

BlP(x)=0 


A2(x)=0 

A2P(x)=0 

B2(x)=-ts2/(e2e0) 

B2P(x)=0. 

As  another  example,  consider  a  one-layer  medium  of  depth  1  m  (i.e.,  zQ=0, 
zt=l)  with 

e(z)=9(z+2)2e0,  0<Z<1 

o(z)=( 10-(z+l)2)»lo”3  mho/m,  0<z<l. 


Then 

?  1 

X=  -f  (s+2)ds=15/2c , 

C  0 

x=x(z)*( z2+4z)/5 
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and 


z=z(x)=-2/$+5x. 


It  follows  that 

Al(x)=2.5/(z+2)2=2.5/(4+5x) 

AlP(x)=-12.5/(4+5x)2 

Bl(x)=-10”3 (5+5x-2/?+5xy/(6ce0{4+5x} ) 

5/2 

BlP(x)=-5* 10  (4+5x-/5+57y/( 6c eQ {4+5x}  ). 

For  more  general  permittivity  profiles,  it  may  be  necessary  to  solve 
equation  (B.l)  numerically  for  x  as  a  function  of  z  using,  say,  Newton's 
method.  More  will  be  said  about  this  in  Section  VI. 


III.  INPUTTING  THE  INCIDENT  HAVE 
The  "time"  variable  t  used  in  HATS  has  been  normalized  by  settinq 

T=t/£ 

and  hence,  t  is  a  dimensionless  quantity.  Thus,  if  the  incident  field  is 
E^(-t),  t>0,  then  the  corresponding  incident  field  UINC  to  be  used  In  the  pro¬ 
gram  is 


UINC(-t)*E1(-t£),  t>0 


or 


UINC(t)*E1 (it) ,  T<0. 
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For  example,  consider  an  incident  sinusoidal  wave  train  of  frequency  f.  Then 


E1  (-t)=sin(2irft) ,  t>0. 

It  is  assumed  that  E1(-t)=0  for  t<0.)  The  incident  field  to  be  used  in  the 
program  is 

UINC(T)=sin(-2nftT),  t <0 • 

This  function  is  to  be  inserted  into  the  subroutine  INCWAV.  Depending  on  the 
complexity  of  the  incident  field,  a  single  function  statement  may  not  suffice 
to  define  the  incident  wave.  This  calls  for  an  obvious  modification  of  the 
subroutine.  See  the  source  listing  in  Section  IX  for  an  example  of  such  a 
modification.  In  that  example,  a  wave  train  of  frequency  f  and  of  finite 
duration  is  generated. 


IV.  OTHER  INPUT  DATA 

In  addition  to  the  function  statements  that  must  be  written  into  the  pro¬ 
gram,  certain  constants  are  to  be  read  in  as  data. 

Line  1:  N,  NL,  IFRQ,  ITM,  ISCAT,  I0PT 

Format:  (615) 

N:  number  of  subintervals  to  be  used  in  the  x  interval  [0,1].  Since 
HATS  solves  partial  differential  equations  and  integral  equations 
numerically,  the  x  interval  [0,1]  is  split  into  N  equal  subinter¬ 
vals  to  establish  a  mesh  of  grid  points.  Thus,  H=l/N  is  the  step 
size  used  in  most  of  the  calculations.  To  reduce  the  effects  of 
roundoff  error,  it  Is  suggested  (but  not  necessary)  that  N  be  a 
power  of  2,  say  N=16,  32,  64,  or  128.  N  must  be  an  even  integer 
not  exceeding  128.  Generally,  the  larger  N  is,  the  more  accurate 
the  program  results,  but  this  Is  at  the  expense  of  increased  com¬ 
putation. 
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NL:  number  of  layers  In  the  medium.  NL  must  be  at  least  1  and  no  greater 
than  5. 

IFRQ:  frequency  of  observation  of  the  internal  field.  In  addition  to  using 
a  "spatial"  (x)  step  size  of  H,  HATS  also  uses  a  "temporal"  (t) 
step  size  of  H.  The  internal  field  will  be  output  at  t  intervals 
of 

At=2H* IFRQ; 

i.e.,  at  time  Intervals  of 

At=2H» IFRQ»t. 

ITM:  maximum  "time"  to  be  considered.  ITM  is  the  maximum  number  of  t 

intervals  of  length  2H  to  be  used  in  the  program;  i.e.,  the 
maximum  time  t  considered  will  be 

t=2H* ITM*i. 

ITM  must  not  exceed  2000  and  should  be  an  Integral  multiple  of 

IFRQ. 

ISCAT:  output  option.  If  ISCAT=1,  the  incident,  reflected,  and  transmitted 
fields  (i.e.,  field  transmitted  out  the  right-hand  edge  of  the 

medium)  will  be  printed  out  for  t=0,  2H,  4H,  •••,  2H»ITM.  If 
ISCAT*0,  this  printout  Is  omitted. 

IOPT:  output  option.  If  I0PT*0,  the  field  within  the  medium  will  be 

printed  out  at  the  (l+N/2)  grid  points,  x-0,  2H?  4H,  •••,  1  for 
values  of  t  given  by 

t=2H» IFRQ,  4H.IFRQ,  6H»IFRQ,  ...»  2H*ITM. 

If  I0PT*»1,  the  field  will  be  printed  out  for  selected  qrld  points 
(to  be  supplied  later)  at  the  above-mentioned  t  values. 


Line  2:  X(l),  X(2) 


,  X  ( 6 ) 


Format: 


Line  3: 

Format: 


(615) 

Assuminq  that  layer  i  extends  from  z=z^ _ ^  to  z=z.j ,  the  Liouville 
transform  (B.l)  implies  that  in  the  x  variable  this  layer  extends 
from  x^_j  to  x^,  where 

zi  i 

x1-(/  {g(s)m0} ds)/i,  i=0,  1,  NL. 

The  quantities  X(i)  are  defined  by 


X(i)=xi_1/H. 


Notice  that  X(1)=0  and  X(NL+1)=N.  The  X ( i )  * s  are  to  be  even 
integers  satisfying 

X(i)-X(i-1)<32,  i=2,3,  •••,  NL+1. 

If  the  medium  has  fewer  than  five  layers,  the  values  assigned  to 
X(NL+2),  •••,  X ( 6)  are  immaterial. 

C(l) ,  C(2) ,  C ( 6 ) 

(6E12.5) * 

If  the  permittivity  in  layer  i  is  given  by  e.(z),  then 

C(1)={ci(zi+)/ef_1(zi-)}i  (B.2) 

where,  as  usual,  the  z^'s  denote  the  boundaries  of  the  layers. 

For  1*1,  (B.2)  is  replaced  by 
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1 


and  for  1=NL+1, 


C(l)-{e1(0+)/c(J}i 


C(NL+l)-{e0/eNL(L-)}*. 


Again,  If  there  are  fewer  than  five  layers,  the  values  assigned  to 
C(NL+2),  •••,  C(6)  are  Immaterial. 

The  following  data  must  also  be  supplied  if  I0PT=1. 

Line  4:  NOBPTS 

Format:  (15) 

If  It  Is  desired  to  determine  the  internal  field  only  at  selected 
spatial  grid  points,  then  NOBPTS  Is  the  number  of  points  in  the  x 
interval  [0,1]  at  which  the  field  should  be  determined. 

Line  5:  IX 
Format:  (15) 

IX  Is  an  Integer  specifying  the  position  of  the  first  grid  point 
at  which  the  field  should  be  determined.  It  Is  defined  by 

IX-x/H 

where  x  Is  the  point  where  the  solution  Is  desired.  IX  must  be  an 
even  integer,  QlIX^N. 

Line  6:  IX 
Format:  (15) 

This  specifies  the  position  of  the  second  point  of  observation. 
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This  same  input  is  continued  until  all  NOBPTS  positions  have  been  speci¬ 


fied. 

V.  OUTPUT 

HATS  automatical ly  prints  out  the  input  data  from  lines  1,  2,  and  3. 

If  ISCAT=1,  the  incident,  reflected,  and  transmitted  fields  are  printed 
next.  The  incident  field  is  printed  first  in  (1X,8E10.3)  format  for  as  many 
lines  as  is  necessary  to  print  the  entire  field.  The  array  being  printed  out 
is 

UI(J),  J=l,  2,  •••,  ITM+1 

where 

MI (J )=UINC(-2( J-l )H) 

=E1 ( —2 ( J— 1 )HA ) 

The  reflected  field  is  printed  next,  using  the  same  format.  The  array 
being  printed  here  is 

UR(J)=Er(2(J-l)Ht),  J-l,  2,  •••,  ITM+1. 

Finally  the  transmitted  field  is  printed  where 

UT(J)=Et(-2(J-l)Ht),  J-l,  2,  •••,  ITM+1. 

The  remaining  output  depends  on  the  value  of  IOPT.  If  I0PT=0,  the  next 
line  of  output  is  an  integer  IT  (in  110  format)  followed  by  an  array  of 
(l+N/2)  numbers  in  (1X,8E10.3)  format.  The  array  being  printed  out  is 

U(J)=u(2(J-l)H,  2* IT*H) . 
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In  other  words,  the  field  throughout  the  width  of  the  medium  is  being  printed 
for  time  t=2*IT*H.  To  relate  this  to  the  physical  field  E,  note  that 


E(z,t)=u(x,i) 
where,  if  x  ._j<x<x,. , 

!+(*JX  exp{-/y  A(s)ds}dy)/(e.(zi_1+)li0)\ 
xi-l  xi-l 


t=x£ 


and 


A(s)=Ai(s)  for  xul<s<xv 

This  printout  is  continued  for  all  values  of  IT  given  by 
ITMFTQ,  2.IFRQ,  •••,  ITM. 

On  the  other  hand,  if  I0PT=1,  then  the  next  line  of  output  is  an  integer 
IX  (in  110  format)  followed  by  an  array  of  l+(lTM-(IX/2))/lFRQ  numbers  in 
(IX, 8310. 3)  format.  The  array  being  printed  is 

U(J)=u(IX*H,  2( IFRO* (J-l )+IX/2 )H) . 

In  other  words,  the  field  at  position  x=IX»H  is  being  printed  at  t  inter¬ 
vals  of  2H» IFRQ  starting  from  t=IX*H.  This  printout  is  continued  for  all 
NOBPTS  values  of  IX  specified  by  the  user  in  the  input  data. 

Finally,  the  message  FINISHED  Is  printed  to  indicate  that  the  program  has 
finished  execution  in  a  normal  fashion. 

In  general,  deducing  any  reasonable  conclusion  by  looking  at  the  printed 
output  will  be  a  formidable  task.  The  user  will  probably  find  it  much  more 
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meaningful  to  plot  the  data,  using,  for  example,  a  CALCOMP  plotting  routine. 
For  this,  a  two-step  process  is  suggested. 


1.  Run  HATS,  havinq  the  output  transferred  to  a  disk  file. 

2.  Run  the  plotting  routine,  readinq  data  from  the  disk  file. 

To  facilitate  step  1,  all  data  is  output  from  HATS  via  the  subroutine 
PRNT.  The  comment  statements  in  the  subroutine  are  self-explanatory,  indicat¬ 
ing  which  blocks  of  code  are  used  to  output  the  various  data.  In  addition, 
the  commented  write  statements  can  be  used  to  send  a  streamlined  form  of  the 
output  to  a  disk  file.  The  user  should  note  the  format  used  so  that  the  data 
is  later  read  properly  from  the  disk  file.  Also  note  that  output  device  8 
needs  to  be  defined  via  JCL. 

An  example  of  a  plotting  routine  has  been  listed  in  Section  X.  The 
actual  calls  made  to  the  graphing  subroutines  will  have  to  be  replaced  with 
calls  utilizing  software  at  USAFSAM. 


VI.  PRECAUTIONS,  TRICKS,  AND  GENERAL  OBSERVATIONS 

As  with  all  numerical  solutions,  HATS  is  not  infallible.  Continuous  but 
wildly  vary inq  conductivity  and  permittivity  profiles  within  a  layer  will  gen¬ 
erally  cause  inaccurate  solutions  to  be  calculated.  Also,  if  the  conductivity 
is  excessively  large,  the  numerical  solution  of  the  partial  differential  equa¬ 
tions  will  aqain  prove  unreliable.  In  this  case  Internal  and  scattered  fields 
will  generally  start  to  grow  with  time,  so  it  will  quickly  become  apparent  to 
the  user  that  the  results  are  in  error. 

One  method  for  trying  to  handle  large  conductivities  or  wildly  varying 
conductivities  and/or  permittivities  is  to  decrease  step  size  H  (i.e., 
increase  N).  Since  this  will  generally  cause  the  restriction 

X  ( 1 ) -X  ( 1  -1  )_<32 

to  be  violated,  the  following  example  shows  how  HATS  can  be  "tricked." 


Consider  a  two-layer  medium  and  suppose  N=32,  with  X(1)=0,  X(2)=12, 
X ( 3 ) =32 .  If  a  smaller  step  size  Is  needed  and  N  should  he  64,  this  will  place 
the  boundaries  of  the  layers  at  X(1)=0,  X(2)=24,  X(3)=64;  consequently 

X(3)-X(2)=40>32. 

To  circumvent  this  problem,  define  the  medium  to  consist  of  three  layers 
instead  of  two,  with  boundaries,  for  example,  at  X(l)=n,  X( 2) =24 ,  X(3)=44, 
X(4)=64.  The  function  statements  for  A3(x),  A3P(x),  B3(x),  R3P(x)  will  be 
exactly  the  same  as  those  for  layer  2.  Further,  C(l)  and  C(2)  will  be  as  in 
the  N=32  case,  and  C(4)  will  have  the  value  that  C(3)  has  in  the  N=32  case. 
Finally,  to  indicate  that  there  Is  indeed  no  discontinuity  at  X(3),  set 
C ( 3 ) = 1 .0. 

If  N  needs  to  be  increased  beyond  128,  consult  the  last  portion  of  this 
Section. 

If  a  large  step  size  (small  N)  is  used,  then  detail  in  the  solution  is 
lost  and  some  accuracy  in  the  solution  values  calculated  will  generally  be 
lost.  On  the  other  hand,  fewer  calculations  will  be  performed  within  the  pro¬ 
gram  and  the  solution  can  be  monitored  over  a  longer  period  of  time  since  the 
maximum  value  of  t  is  ITM*2H. 


Extending  the  Maximum  Time 

If  necessary,  values  of  ITM  greater  than  2000  can  be  considered  rela¬ 
tively  easily,  at  the  expense  of  increased  storage.  It  suffices  to  increase 
the  dimension  of  the  arrays 

V,  BUFF,  UI,  U,  UR 

which  appear  in  COMMON  statements.  Let  ITM1=ITM+1.  The  new  dimensions  for 
the  above-mentioned  arrays  are  to  be 

V(5,ITM1,2) ,  BUFF(200+ITM1 ) ,  UI(ITMl),  U(ITM1),  UR(ITMl). 


40 


Incidentally,  if  there  are  fewer  than  five  layers,  storage  can  be  reduced  by 
changing  the  5  in  V(f>,ITMl,2)  to  the  actual  number  of  layers.  (In  subroutine 
FIELD  there  is  an  array  UT(2001)  that  contains  the  transmitted  field.  The 
dimension  of  this  array  need  not  be  increased  since  UT  is  eouivalenced  to 
BUFF). 


Increasing  the  Number  of  Intervals  A 1 1  owed 


If  a  value  of  NL  greater  than  5  is  desired,  it  is  necessary  to 

1.  Change  the  dimension  of  all  arrays  in  COMMON  statements  from  5  to  the 

new  value  of  NL  (i.e.,  wherever  a  5  appears,  it  should  be  replaced 
by  NL). 

2.  Change  the  appropriate  READ  statements  in  the  main  program  and  WRITE 

statements  in  subroutine  PRNT. 


Decreasing  the  Step  Size 

If  a  value  of  N  greater  than  128  must  be  considered,  it  is  necessary  to 

1.  Change  the  dimension  of  all  arrays  appearing  in  COMMON  statements 

from  129  to  N+l. 

2.  Change  the  dimension  of  BUFF  to  N+1TM+3. 

3.  Change  the  EQUIVALENCE  statement  at  the  beginning  of  subroutine  FIELD 

to 

EQUIVALENCE  (BUFF(N+2),  UT(1)). 
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Numerical  Solution  of  Equation  (B.l) 


For  some  permittivity  profiles  the  Llouvllle  transform  (B.l)  will  need  to 
be  solved  numerically  for  z  as  a  function  of  x.  This  requires  two  new  arrays 
to  be  generated,  Z(I)  and  ZH(I),  1=1,  2,  ••*,  N,  where 

Z(I)  i 

(I-l)H=(Jo  {e(s)u0}  ds)/l 
ZH(I)  j 

(I-.5)H=(/o  { e ( s ) M0 }  ds )/i 

and  Z(N+1)=L.  The  above  equations  are  readily  solved  via  Newton's  method. 

Since  the  functions  AI,  AIP,  BI,  BIP  In  subroutine  EVAL  cannot  be  ex¬ 
pressed  In  terms  of  the  x  variable,  they  must  now  be  written  as  functions  of 
z. 


j 

AI(Z)='  dl{e  i(z)mo>  ’  zi_i<z<*i 
AIP(z)=(±JU(z))/(£i) 

OZ  QZ 

BI(z)=-JlOj(z)/eI(z),  zN1<z<Zj 
BIP(z)=(zl8I(z))/(cI(z)M0}i. 

Finally,  In  subroutine  EVAL  the  statements 

XR=(I-1)*H 

XR=(I-.5)*H 

should  be  replaced  with 

XR«Z(I) 

XR=ZH(1) 
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respectively,  and  the  remaining  eight  statements  of  the  form 


r 


XR=(X(J)+I-1)*H 

XR=(X(J)+I-.5)*H 

should  be  replaced  with 

XR=Z(NPTL(J-1)+I-1) 
XR=ZH(NPTL(J-1)+I-1)  • 


VII.  SOURCE  LISTING  OF  HATS 


..........MAIN  PROGRAM 

INTEGER  X 

CC  XMCN  /  BLOCK  1/11.  NL  «  NTGP  , NPTL (  5)  .  I  FRO.  ITM.  IT.  I  SCAT.  I  OPT 
COMMON  /3LOC*2/A<  331.8(33)  .AP( 33 » , BP <  33) . DT { 3 3  I . AH ( 3  3 )  ,BH< 33) . 

1  APH<33) .3PM (33) »OTH(23 )  »  A  J ( 33 ) . 31 ( 33> . 3T I ( 33 ) . 

2  OP  (  33  )  .  cB  I  (  33  )  .  AH  (  5  )  .  AO  (  5  >  .  BR  (  5  )  .  BL  (  5  ) 

COMMON  /dLOCK3/X(6).C(6).OPI  (  3  3  )  ,  DPB I  (  33  >  .  E  A  MS  ( 5 . 33 )  ,£A»B(S.33) 
>  UC.LDI  33  )  ,  ONEX  (33  )  *XP(  S.  129  )  ,  XM<  5*1  29 1 .  YP<  5.12V). 

2  YM(5. 129) .XPC(5) . XMC ( 5 ) . YPC( 5) . YMC( 5) . APB 1 ( 5) . 

3  AMBl(S) ,P( 5.2.33.33) 

COMMON  /3L0CK*/UV(  S,16).J«r(5.16).SV(S.le)  ,  SW(  5,16)  ,V<  S.20ul  ,2)  . 
1  W<  5.129.2). V5( 129.2 ).WS( 129,2) 

COMMON  /BLOCKS/dUFF ( 2201 >,UI (2001 ) *01200 l ) .OR (2001 ) 

READ  (5.2000)  N.  NL,  IFRu.  ITM,  ISCAT.  I  OPT .  (X(I),I  *  l.o) 

READ  (5.3000)  ( C( I ) , I  =  1.6) 

NT  GP  *  n  4  1 

M  *  1 .0  /  N 
CALL  PRNT( 1 ) 

CALL  KERNEL 
CALL  DELTA 
CALL  INC WAV 
CALL  FIELD 
CALL  PRNT16) 

STOP 

2000  FORMAT  (615./. 615) 

3000  FORMAT  ( 6E12.5 ) 

END 


SUBROUTINE  E VAL(LAYER . ICK) 

.....INSERT  FUNCTION  STATEMENTS  FOR  EACH  LAYER 

A1 (X)  *  0.0  *  0.  •  X 

31 ( X)  *  0.0  *  0.  *  X 

A 1  P(  X )  •=  0.0  ♦  0.  »  X 

31P(X)  =  0.0  *  J .  »  X 

A2 ( X )  *  0.0  ♦  0.  *  X 

62 ( X )  *  0.0  4  0.  *  X 

A2P( X)  *  0.0  4  0.  *  X 

32P{ X )  *  0.0  4  0.  *  X 

A3 ( X )  a  0.0  4  0.  *  X 

33 (X)  a  0.0  4  0.  •  X 

A3PCX)  s  0.0  ♦  0.  *  X 

33P(X)  9  o.o  4  0.  *  X 

AA(X)  a  0.0  4  0.  •  X 

3A(X)  a  0.0  «•  0.  *  X 

AAP(X)  a  0.0  4  0.  *  X 

B*P(X)  a  0.0  ♦  0.  *  X 

A5CX)  a  0.0  ♦  0.  *  X 

35(X)  a  0.0  4  0.  *  X 

ASP(X)  a  0.0  4  0.  *  X 

B5P( X )  a  0.0  4  0.  •  X 


INTEGER  X 

COMMON  / BLOCK 1 /H» NL. NTGP. NPTL ( 5) .IFMO. ITM • IT . I SCAT • I OPT 
CO  MMON  / BL0CK2/ A ( 33 ) . B ( 33 ) • AP<  33 ) t  BP <  3  3 ) , OT  <  33  > • AM  < 33 ) . BH( 331. 

1  APH(33>.BPH(33).OTH(23).At (33) .SI ( 33) . O T I ( 33 ) . 

2  OP(33).Eai(33).AR(3)4AL(6)*SR(6).BL(S) 


89  • 
90. 
91* 
92. 
93* 
94.. 

95. 

96. 

97. 

98. 

99. 

too. 

101. 

102. 

103. 

104. 

105. 

106. 
107. 

ioa. 

109. 
1  10. 
ill. 
1  12. 
1  13. 
I  14. 
1  15. 
1  16. 
I  17. 
118. 

119. 

120. 
121  . 
122. 

123. 

124. 

125. 

126. 

127. 

128. 

129. 

130. 

131. 

132. 

133. 

134. 

135. 

136. 

137. 
133. 

139. 

140. 

141. 

142. 

143. 

144. 


OT  Cl)  *  (Bd)  M2*  Adi  •*  21  /  4.0 
X*  *  (X(4l  41*  .3)  4  H 
AH  C  1  )  ■  A4(XR) 

BH(l)  *  B4CXAJ 
APH(Z)  *  A4P<XR) 

BPHO)  *  B4P(XR) 

OTHd)  *  <BH(1)  mm  2  -  AH(  I )  M  2)  /  4.0 
55  CONTINUE 
GO  TO  67 

60  00  65  1  *  1.  NP1 

XP  *  <X<5>  ♦  I  -  I .  I  OH 
Ad  )  *  A 3(XR) 

B< 1 )  «  85C  XR1 
A Pd)  m  AS PCXR1 
BP  (I )  *  8SPCXR) 

OT  (II  »  (0d)  *02-  Adi  *•  2  i  /  4.0 
XR  *  ( X( 5 I  ♦  X  -  .5)  •  H 
AH  d )  *  ASCXR! 

BHd)  *  85 C  XR  ) 

APH(l)  «  ASP(XR) 

BPH( 1 )  *  85PCXR) 

OTH(I)  *  (BHd)  *o  2  -  AH  CXI  OO  2)  /  4.0 

65  CONTINUE 
67  CONTINUE 

AI  CNPi  )  *  0.0 

ai (npi j  *  o.o 

0TKNP1)  »  0.0 
EB I <  NP 1)  *  1.0 
N  »  NP 1  -  1 

DC  68  K  =  1 t  N 

J  =  NPI  -  K 
JP  *  J  ♦  1 

AI(J)  a  A 1  (  JP )  «■  M9  m  (A  (J)  4  4.  4  AH(J)  4  A(JP)| 

ai(J)  «  3KJP)  4  H6  *  <a(J>  4  4.  O  BHtJl  4  3(jP)> 

DTI(J)  s  OTI(JP)  ♦  HS  *  1DT1J1  ♦  4.  *  OTH(J)  ♦  OT  t  JP )  ) 

EEI(J)  a  €XP (  BI(J)  ) 

66  CONTINUE 
GO  TO  90 

70  OC  80  1*1.  NPI 

Od  )  *  —  3(1) 

8Pt I i  a  -  8PI 1 } 

BH( 1 )  *  *  3Ht I ) 

BPHd)  a  -  8PH(Ii 
Bid)  «  *  Bid  ) 
esid  )  a  1.0  /  £31  <  1  i 
•0  CONTINUE 
90  CONTINUE 

AR (LAYER )  a  At NPI ) 

0fi (LAYER)  *  8CNPI) 

AL( LAYER)  a  Ad  ) 

BL (LAYER)  •  8(1) 

00  100  I  «  li  NPI 

0P(1)  *  0T(I)  *  (APd)  *  BPdl)  /  2. 

100  CONTINUE 
RETURN 
END 
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1*  SUBROUTINE  INC*AV 

C« .........  INSERT  ^-UNCTION  STATEMENT  for  INCIDENT  FIELD 

3  •  J I NC (X)  s  0  •  0  ♦  0 •  *  X 

4.  . . 

5*  COMMON  /dLOCKi/H.NL.NTGP .NPTLI  5> .IFWU. I TR , I T. I  SCAT. I  OPT 

°*  COMMON  /  BLOCKS /BUFF  (2201 )»UI 12001 ).U(200l  )  .  UR  (  20  01) 

T.  H2=H*2. 

9  •  ITM1  =  ITM  ♦  1 

9.  DO  10  J  =  1 «  ITM I 

10.  X  =  —  (  J  -  1  )  *  H2 

11 •  UI I J)  =  UINC(X) 

12.  10  CONTINUE 

13.  RETURN 

14.  END 

1.  SUBROUTINE  PRNTIK) 

2.  INTEGER  X 

3.  COMMON  /BLDCK1/H. NL .NTG= .NPTL1  5)  . IFRO. I TA . I T, I  SCAT .IOPT 

A.  COMMON  / BL0CK2/AI 33 ) .3(33 ) • API  33 ) • BP (331* DTI  331*  AHt  33 ) »SH( 33)  . 

5.  1  APHI 33)  .BPH 133)  .DTHI 33) . All  33 ) ,81 ( 33 ) .OTI  < 33 j , 

6.  2  l)PI33)  .EBK  33).  ARIS)  .  ALIS)  .BR(  5).BLI5) 

7.  COMMON  /BLJCK3/XI 6) .Cl  6) . DP I  I  33 )* OPS I  I  33  > . EAMB 15.33) .EAPBI 5.33) . 

B.  1  QOLDI 33 ) ,0NE«I33) , XP I  5 . 1 29 ) . XM I  5. 1 29 ) .VPI5.129) . 

9.  2  VMI 5. 129) .XPCI5) . XMCC 5) . YPCI 5) . YMCI 5) . APBI l 5) . 

10.  3  AMBUS)  .PIS. 2.33. 33) 

1 1 .  COMMON  /3L0CK4/JVI 5. 16) • JWI S.l 6) .SVC S • 16) .S*( 5* 1 6 ) .V  I  5.2001 .2) . 

12.  1  w( 5.129.2) . VS< 1 29.2 ) .MSI  129*2) 

13.  COMMON  /BL0CK5/BUFFI 22  01  )  .UII2001  )  .1)1200 1  )  , UR  1  20  0  1  > 

14.  IF  IK  -  2)  1.  2.  SO 

15.  I  CONTINUE 

IS.  C.. ..... ...PRINT  OUT  INPUT  DATA  FROM  LINES  1.  2.  3 

17.  N  =  NT GP  “1 

IB.  NLL  =  NL  ♦  1 

19.  WRITE  16.2000)  N.  NL.  IFRQ.  ITM,  ISCAT.  I  OPT ,  IXII).  I  *  1.  NLL) 

20.  WRITE  16.3000)  <C(I).I  =  1.  NLL) 

21.  C  WRITE  IB. 2000)  N.  NL.  IFRO.  ITM 

22.  RETURN 

23.  2  CONTINUE 

24.  RETURN 

25.  50  CONTINUE 

26.  IF  IK  -  4)  3.  4.  60 

27.  3  CONTINUE 

23.  C. ....... ..PRINT  INCIDENT,  REFLECTED  ANO  TRANSMITTED  FIELOS  IF  ISCAT  =  1 

29.  IF  I ISCAT  .EO.  0>  RETURN 

30.  ITM!  *  ITM  ♦  1 

31.  WRITE  16.5000) 

32.  WRITE  16.1000)  (UI(I).I  «  1.  ITMt  ) 

33.  '  WRITE  1 6. SI 00 ) 

34.  WRITE  16.1000)  I URt 1 ) . 1  *  1.  ITM1) 

35.  WRITE  16.5200) 

36.  WRITE  (6.1000)  (V(NL.l.l).I  »  1.  ITMt) 

37.  C  WRITE  (8.1500)  (V(NL.l.l).I  *  1.  ITM1  ) 

38.  C  WRITE  (8.1500)  (UR(I).I  «  1.  ITM I ) 

39.  RETURN 

40.  4  CONTINUE 

41.  C... ...... .PRINT  FIELD  AT  ALL  GRID  POINTS  INSIDE  MEDIUM  AT  TIME  TAU  *  1  * 

42.  WRITE  (6,7000)  IT 

43.  WRITE  (6.1000)  (U(I).I  *  1.  NTGP.  2) 
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44  * 

45  * 
46. 
4  7  * 
4  6  * 

. 


-  . 

5  1  . 


C 

c 


60 


C. 

C. 


•  RITE 
WRITE 
RETURN 
CONTINUE 
IF  (K  -  6 i 
CONTINUE 
,  •• •••print 
. . . 


(3*6000) 
(  3*1500) 


IT 

(U(I >*I 


»  i*  NT  GP •  2) 


f i Eld  AT 

»  XT  •  H 


POSITION  X  * 
.  (IT  9  IFRQ 


IT  *  H 

*  2)  *  H, 


FOR  TIMES 


52. 

IMAX  * 

i  -#■  citm  -  (IT  /  2)  J  /  ing 

53  * 

•rite 

(o.aooo)  it 

54  . 

•RITE 

(6.1000)  (U(l).l  *  1.  IMAX) 

55  * 

C 

•  RITE 

(a. 6000)  IT 

66. 

5  7  . 

C 

•  RITE 
RETURN 

(S.1500)  (U(I).l  *  1.  IMAX) 

5o  « 

b 

continue 

3  ^  A 

6  3  . 

•  RITE 

RETURN 

(0*40GO) 

6  l  . 

It  00 

format 

( ix.a&io.3) 

52 

15C0 

FORMAT 

( 3E10.3) 

t  5  - 

2C  0  0 

format 

(615./. 615) 

0*  * 

30L  0 

format 

(6E12.5) 

4COO 

FORMAT 

( ///. 1 X. SHF  INI SHE01 

6  * 

50  0  0 

FORMAT 

(///.1X.14HINCIOENT  FIELD*/) 

67* 

51C0 

FORMAT 

(///* IX* 15HREFLECTED  FIELD*/) 

65  « 

5200 

format 

(///. IX.17HTRANSMITTED  FIELD./) 

o  c>  * 

5000 

FORMAT 

(110) 

7  0  . 

7000 

FORMAT 

(///. 1X*29HINTE«NAL  FIELD  AT  T I Mt  TAU  *  • 

t  l  4 

7  2  » 

acoo 

FORMAT 

END 

(///♦ 1X.21HF1EL0  AT  POSITION  X  A.110.4H  A 

•  2  *  H./ 


«•/  ) 


1  . 


6  % 

7  . 

5. 

<5. 

10. 

1  1  . 
12. 
13* 
14* 
16. 
10. 
17. 
16. 
1  ** 
20. 
21  - 
22. 

23. 

24. 
2S- 
2c. 
27. 
26. 


SUBRCUTINE  kernel 

. . this  subroutine  solves  SYSTEMS  of  PARTIAL  DIFFERENTIAL 

c . * . . . . EDJATICNS  TC  OBTAIN  RIEMANN  FUNCTIONS  FOR  EACH  LAYEh 

. . OF  THE  MEDIUM 

INTEGER  X 

COMMON  /BLOCK 1/H. NL • NTGP*  NPTL ( 5) « IFRQ* ITM*IT*ISCAT*I0PT 
COMMON  /3L0CK2/AC 33) * 0 ( 33 ) * AP( 33  > • 0*t  33) * OT 1 33 ) *  AH( 33 ) • BH( 33 ) • 

1  APH( 33 ) ♦ 0PH ( 33 ) *DTH( 33 ) *AI (33 ) *01 ( 33) • DTI (33 ) . 

2  OP (33  I • £01(33). ARC 5 )*AL(5> *0R(5)*0L(5> 

COMMON  /0LOCK.3/X  (  6)«C(6l*0PI(33)  •  ORB  I(  33  )  *EAM04  5*  33)  .EAPBf  5*33)  » 

1  QOLOC33) t ONEWC  33) *XRC5? 129) *XM( 5* 129) » YP( 5* 129) * 

2  TMC  5* 129 ) ♦ XPC C  5 ) * XMC (  5 )  * VPC( 5) • VMC( 5) *AP0Ii 5) • 

3  A MB 1(5) *P( 5*2 *33*33) 

DIMENSION  W  (33 ) 

C. ...... ...NL  a  NUMBER  OF  LAYERS 

C. ..*•** •••NPTL( I )  *  NUMBER  OF  GRID  POINTS  IN  LAYER  I 


00  1 0  I  «  1  *  NL 

NPTL<  I)  *  XC1+U  •  X(  I  )  9  1 
10  CONTINUE 
H4  *  H/4 
HS  *  H/6 

00  110  I  •  1  *  NL 

FI  m  •  CCII91)  ♦  I • I  /  Z*0 

F2  *  (C( 1911  9  1#)/  6#0 

F3  *  (  C  (  I  9 1  )  -  1  •  )  /  2*0 

F4  a  F3  /  4.0 

F5  «  F3  /  2*0 

NP1  *  NP  TL ( I ) 
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29.  N  »  HP  1  -  I 

30*  DO  100  ICK  «  1 .2 

31.  CALL  £  VAL (  1  •  ICK) 

32*  F6  ■  (ARCI)  -  8R<  X ) )  /  4*0 

33.  01  a  F3  *  C  AR (11  ♦  8R(I)>  ✓  2.0 

3* •  02  «  •  F2  l  F6  *  4*0 

35.  03  *  F2  4  OP(NPl) 

36.  0PKNP1)  s  0*0 

37.  DPQKNPl)  *  0.0 

38.  AMQ  *  4  AR (XI  -  BR 4  I  j J  /  2.0 

39.  DO  20  K  at  1 .  N 

40*  J  k  NP1  *  K 

*!•  JP  *  J  ♦  l 

42*  OP  U  J )  *  DTI(J)  ♦  (ACJ)  -  B(  J }  )  /  2.0  -  AMB 

43*  0P31  ( J  )  *  0P8KJP)  ♦  HS  •  (  B(  J  )  *  4  DT  (  J  >  -  ( AP{  J  >  -  BP(JJ)  /  2.0> 

44*  i  ♦  4.  *  BH(J)  *  CDTH(J)  -  (APH(J>  -  dPH<J))  /  2-0) 

2  ♦  9(JP)  4  (OT(JP)  -  (AP(JP)  -  BP(JP))  /  2.0)) 

4°*  20  CONTINUE 

47*  C.  •••••••* .COMPUTE  BOUNDARY  VALUES 

4S.  DC  25  K  a  1,  NPi 

49*  P(X.ICK.K.l)  »  Fi  *  (OTI(K)  ( A  ( K )  -  BCK)>  /  2.0)  4  01 

-0.  J  *  NP1  ♦  1  -  K 

5l*  P(X.lCK.NPl.K)  *  EBI(J)  A  (F3  *  (OTI(J)  ♦  IA(J)  4  B(J))/2.0)  ♦  OZi 

52.  25  CONTINUE 

53.  .  C. •••••• ...GENERATE  SOLUTION  OF  PQE 

54.  00  40  NCT  a  2.  N 

55.  LN  =  N  ♦  2  -  NCT 

56.  KM  A  X  =  NCT  -  1 

5T.  ON  s  DP{ LNJ 

55.  DNM  =  DP ( LN- | ) 

59.  EN  =  EBI (LN ) 

60.  OL LO ( 1 )  a  F2  •  ON/ EN 

61  .  5NM  s  EDI  (LN-1  ) 

62.  JCLO(NCT)  =  J3  +  <F4  a  OP  I(LN)  -  F6)*0PULN)  ♦  F5  •  DP  31  ( LN ) 

63.  EQ  =  ENM/EN 

C4.  H4  s  H / 4  . 

65.  HE  =  H/EO 

66.  DE  —  ON/cN 

67.  HDE  a  H4  *  ON  »  EG 

65.  DET  *1.+H4*H*  DNM 

69.  DO  30  J  a  1,  KMA  X 

70.  P  *  LN  ♦  J 

71.  Mi  a  M  -  i 

72.  J1  *  J  4  l 

73.  P< I .ICK. Ml . J1 )  *  (P< I .ICK.M. J1 )  -  H  4  < NOE  *  P ( I t I CK • M 1 • J )  ♦ 

74.  1  EN  4  OOLD(Jl)  ♦  ENM  4  QQLO( J1 ) )/OET 

75.  QNEM(Jl)  a  (QOLD(J)  +  H4  4  ( DE  4  P ( I • I CK. Ml • J )  ♦ 

76.  1  ONM  4  ( P( I  *  1 CK* M. 01  ) /ENM  -  HE  4Q0L0C J l  )  ) )  ) /DET 

77.  30  CONTINUE 

73.  DO  35  J  *  2* NCT 

79.  QOLO(J)  a  QNEW4JJ 

60.  35  CONTINUE 

61.  40  CONTINUE 

82.  QNEtf(l)  a  F2  *  DNM/ENM 

63.  ONE*  ( NP 1  )  «  03  ♦  C  F4  4  OPlll*  -  F6)  4  OPICU  ♦  F5  4  DPBI  (  1  ) 

64.  . . COMPUTE  INTEGRAL  FOR  SCATTERING  KERNEL 

65.  F7  a  EBI  41)  4  2. 

86.  F8  «  (AL(tl  -  BUD)  /  2-0 

67.  W<  X)  «  0. 
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80o 

99. 

90  * 

91  « 
9?  o 
<53. 
94  s* 
95o 
96  o 
97. 
98  . 
99. 

I  00  . 
10  1. 
102* 
1  03  o 
104* 
1  C  5  . 
1  06  « 
1  07. 

1  03* 
10  5. 
1  10  . 
I  l  l  * 
112* 
113* 

1  U, 

:  is* 

1  1  O  a 

117. 
115* 
119. 
1  20  * 
121* 
122* 
1  23. 

124. 

125. 


Ml  *  F7  MQNEM(l)  -  F8  *  P(I.ICK.l.l) 

DO  50  K  *  2.  NP I 

M2  *  F7  *QNEi( K )  -  F8  •  P( I • ICK.K.K) 

«(K)  *  M(K-l)  ♦  H  *  (Ml  ♦  M2) 

Ml  *  M2 
50  CONTINUE 

IF  ( 1 CK# £Q«2 )  GO  TO  70 
F9  *  EBI ( 1 i 

FI  1  *  F8  *  (CC 1  +  i )  ♦  !•) 

DO  60  K  *  1.  NP1 

XP(I«K)  a  <  -F1I  ♦  M  <  K ) )  /  EBI(l) 

XM(I.K)  *  FI  1  ♦  2.  *  PCJ«1*K*K)  -  WCK) 

60  CONTINUE 

XMC(I)  «  XM(  1 « NP 1  )  -  F3  *  ( AL(  I )  ♦  BL( I ) )  •  E91 ( 1  ) 
XPC(I)  *  -  XMC<I)  /  ESK  1) 

GO  TO  90 

70  G1  *  F3  *  (  AL.  (  1  )  ♦  3L(  X  )  i  *  EB  I  (  1  ) 

G2  *  -FI  *  (  AL  11)  -  eL(I)> 

DO  80  KT  >  U  N?1 
K  *  NP 1  —  KT  ♦  1  • 

VP(I.KT)  a  Gl  -  2.M  P ( I  *  2  • NP I  * NP l )  ♦  Z.m  P(I,2*K*K)  ♦ 
YM(1,KTJ  a  {  — G1  ♦  2.*  P< I*2*NPl .NPl )  -  tf(NPl)  ♦  M ( «0  ) 
80  CONTINUE 


rPCIIJ  a  51  -  2.*  P( 1 .2.NP1.NP1 ) 
YMC1 I )  *  -  YPC(I)  /  EB 1(1) 

G2 

♦ 

M ( NP 1 ) 

s;0  CONTINUE 
*00  CONTINUE 

DO  10S  J  »  1 •  NP1 

EAMB(I.J)  a  E  XP(  (AKJ)  *  BI(J)> 

/ 

2  • 

) 

/  2# 

EAPB(I.J)  a  £XP(  (  A 1  (  J  )  -  BKJ)) 

/ 

2. 

) 

/  2. 

105  CONTINUE 

3R (  I  )  =  -  BR(  I  ) 

dL (  I  )  a  -  BL (  I  ) 

AP9IU)  a  2 •  *  EAP5<  I  •  1  ) 
4M5I(i;  a  2.  *  EAMQ(l.l) 
i 10  CONTINUE 
RETURN 
END 


•  {  NP  I  i-  W  (  K.  ) 
/  Edl(l) 
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1.  SUBKC  U  T 1 NE  DELTA 

2.  c . ••••this  subroutine  combines  the  riemann  functions  from 

3.  C  „  .  . . * . SUBROUTINE  KERNEL  TQ  FORM  THE  SCATTERING  KERNELS 

A  •  INTEGER  X 

5*  REAL*6  3 JM l •  SUM2 »  SUM 3 • SUM A 

6*  COMMON  / BLOCK l /H,NL. NTGP • NPTLC  5  I  . I FRO*  XTM. IT*  I  SCAT  *  I  OP T 

7*  COMMON  /BL0CK3/X ( 6 ) .C ( 6 ) .DP I { 33 > » DPB I ( 33 > .EAMBC5.33) • E  APB  ( 5,33 ) • 

3*  1  GQLDC 33 > • QNEWC33) . XPC S, 129 > ,XM( 5. 129)  . YPCS.129)  • 

9,  2  Y  M  ( 5* 129 ) ,XPC(5)« XMC( 5) .YPCC  51 .YMC( 5) • APB  I <  5) • 

13*  3  AMBIC 5) *P< 5.2.33.33) 

11  ,  COMMON  / BL0CK4/ J V (5«16),Jtf(5*l6).5V(5. 16 ) »SbC 5. 1 6 1  * V ( 5 • 200  I • 2) • 

12.  1  W( 5.129.2) . VS< 1 29.2 ) • ft 5 ( 129.2) 

13.  N  a  NL 

14.  NT GP  *  0 

15.  00  1 0  I  =  l  •  NL 

16.  NT  GP  *  NT  GP  ♦  NPTL(I) 

17.  10  CONTINUE 

Id.  NT  GP  *  NTGP  -  NL  ♦  1 

19.  NT  GP 1  *  NTGP  -  1 

23.  DO  30  I  *  l,  nl 

21.  NP2  *  NPTL(I)  ♦  1 

22.  DO  20  K  a  NP2.  NTGP 

23.  XP( I .K)  *  XPCl I) 

24.  XM(l.K)  *  XMC ( I 1 

25.  YP ( I • K )  =  YPCCI) 

26.  YM(I.K)  *  YMC(I) 

27.  20  CONTINUE 

23 •  30  CONTINUE 

29.  C . CALCULATE  LOCATION  AND  STRENGTH  OF  SINGULARITIES 

30.  NM 1  =  N  —  1 

31.  JV ( N . L )  *  0 

32.  JW(N.l)  *24  X ( N  ♦  1) 

33.  SVCN.ll  *  AM9IIN)  4  CCCN  ♦  11  ♦  1.)  /  2. 

34.  SVCN.ll  *  -  APBI(N)  4  CCCN  ♦  l)  -  1.)  /  2. 

35.  IF  ( N«  EG. 1 )  GO  TO  70 

36.  00  60  IC  *  i.  NM i 

37.  NCL  *  N  -  IC 

38.  NCK  =  NCL  ♦  1 

39.  NJ  *  2  4*  ( N  —  NCL) 

40.  NJ2  =  NJ  /  2 

41.  00  40  I  *  l.  hJ2 

42.  JV(NCL.I)  *  JVC  NCK  .11 

43.  J*(NCL.I)  *  Jm(NCK.I) 

44.  SVCNCL.I)  *  AM3KNCL)  4  ( C C NCK )  4  1.)  4  SVCNCK.il  /  2. 

45.  SW(NCL.I)  *  AP6MNCL)  4  (CCNCKl  ♦  1.)  4  StfCNCK.il  /  2. 

46.  40  CONTINUE 

47.  NJ 3  *  N J2  ♦  1 

48.  00  50  I  *  NJ3.  NJ 

49.  JVC NCL* I  )  *  2  *  X(NCK)  -  JtfCNCK.I  -  N  J  2 ) 

50.  JtfCN/CL.I)  *  2  4  XCNCK1  -  JVCNCK.I  -  NJ21 

51.  SVCNCL.I)  *  -  AMBICNCL)  4  (CCNCKl  -  1.)  4  S«CNCK,I  -  NJ2)  /  2. 

52.  Sir  ( NCL  «  I  )  *  -  AP3ICNCL)  4  CCCNCK)  -  1.1  4  SVC  NCK,  I  -  N  J2  )  /  2. 

53.  50  CONTINUE 

54.  60  CONTINUE 

55.  70  CONTINUE 

3=>.  AP  *  APBICN)  /  4.0 

57.  AM  *  AMBICN)  /  4.0 

56.  DC  80  K  *  1,  NTGP 

59.  VCN.K.l)  *  AM  *  XMCN.K) 

60.  VCN.K.2)  *  VCN.K.l) 
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61.  W(N.K.l)  =  AP  *  XP(N|iO 

62.  W(N.K.2)  *  »(N.K.l) 

£3.  80  CONTINUE 

£*.  V  C  N.  I  *  1  )  *  0. 

65.  Ml N. 1  .  1  )  *  0. 

6o.  NP1  *  NPTL(N) 

67.  VtN.NPl.2)  *  AM  *  XMC(N) 

68.  «(N.NP1.2)  a  AP  a  XPC(N) 

69.  IF  IN. £0.0  CO  TO  2A0 

70.  £>D  23*  1C  *  1.  NMl 

71.  I  ■  N  -  1C 

72.  C....... CALCULATE  INTEGRALS  FOR  LAYERS  N- 1 .N-2, . . .2 *  1 .  IGNORING  SINGULAR  IT 

73.  NP 1  a  NPTLI  I ) 

74.  N1  *  NP  1  -  1 

75.  AN  *  ANSI! I  )  /  *. 

76.  AP  *  APB  I (  II/*. 

77.  IP  »  I  ♦  1 

78.  CM  ■  2.  *  (C(IP)  -  1.1 

79.  CP  •  2.  *  (C(IP>  ♦  1.) 

eo.  VI 1.1.1)  »  0. 

81.  V(I.l.Z)  *  AM  *  (CP  *  VHP. 1.2)  —  CM  *  V(IP*I.2)> 

32.  KCI.l.l)  *  0. 

83.  »( 1.1.2)  *  0. 

8*.  00  lie  K.  a  1.  Nl 

85.  SUM1  a  0.00 

86.  SUM2  •  0.00 

37.  DL  *0  J  ■  l.  K 

es.  IK  a  K  -  J  ♦  2 

89.  JP  a  J  ♦  1 

90.  IL  »  IK  -  1 

91.  SUM1  a  SUM1  *  XM(I.J)  *  VC  IP.  IK.  1)  ♦  XMC1.JP)  •  VC1P.IL.2) 

92.  1  ♦  VMC1.J)  «  W(IP.I«C.l)  ♦  YM(I.JP)  *  pCIP.IL. 2) 

93.  SUM2  a  SUM 2  *  XP(l.J)  •  VC  IP. IK. 1)  ♦  XPCI.JP)  *  VCIP.IL.2) 

9*.  I  ♦  VPCI.J)  *  M<  IP.  IK.  1 )  VPCI.JP)  *  MCIP.IL. 2) 

95.  90  CONTINUE 

96.  KP  •  K  ♦  1 

97.  DO  100  M  a  1.  2 

98.  V(I.KP.M)  a  AM  *  C  CP  *  VCIP.KP.M)  —  CM  a  WC1P.K*».M)  ♦  H  *  SUM1  ) 

99.  WCI.KP.M)  a  AP  *  H  *  SUM2 

100.  100  CONTINUE 

101.  110  CONTINUE 

102.  MCI .NP 1.2)  a  MC1.NP1.2)  ♦  AP  *  (CP  *  MCIP.1.2)  -  CM  *  V(IP.1.2x) 

103.  SUMS  a  0.00 

10*.  SUM*  m  0*00 

105.  00  1*0  K  a  NP 1  *  NTGP1 

106.  SUM1  a  0*00 

107.  SUM2  a  0.00 

108.  00  120  J  »  1.  Nl 

) 09.  IK  «  K  -  J  ♦  2 

l  10.  JP  a  J  *  1 

Ul.  IL  a  IK  -  1 

112.  SUM1  a  SUM 1  *  XMCI.J)  •  VCIP.IK.l)  *  XMCI.JP)  a  VCIP.lL.2) 

113.  .  1  *  VMCI.J)  •  VCIP.IK.l)  *  TM( I .  JP )  •  MCIP.Il.2) 

II*.  SUM2  a  SUM2  *  XPCI.J)  *  VCIP.IK.l)  *  XPCI.JP)  *  VCIP.IL.2) 

I  (5.  1  *  YP(I.J)  •  VCIP.IK.l)  *  VPCI.JP)  a  NCIP.IL.2) 

» 16.  120  CONI  1 MUE 

|17.  So  M3  a  SUM3  ♦  VdP.K  -  NPl  ♦  2.1)  ♦  VCIP.K  -  NPl  ♦  1,2) 

)l3.  SUM*  a  SUM*  *  M(IP.K  -  NPl  *  2.1)  *  NCIP.K  -  NPl  *  1.2) 

119.  KP  a  K  *  1 

120.  KM  a  K  -  NPl  +  2 
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I  z  i  .  DO  i  3  0  M  =  t  *  2 

j  22.  VU.KP.M)  s  AM  •  (CP  *  V(IP.KP.M)  -  CM  *  W(IP.KP.M) 

/  2  3  •  1  4  H  *  (  SUM  1  4  XMC(l)  *  SOM3  4  YMC(I)  *  SUM4J) 

/2*.  4CI.KP.M)  =  AP  *  (CP  *  M(IP.KM.M)  -  CM  •  V (  I P * KM  * M ) 

,25.  1  4  H  *  CSUM2  4  XPC(Ii  *  SUM3  4  YPCII)  *  SUM*)) 

)  26 •  130  CONTINUE 

127.  1*0  CONTINUE 

/  23  *  C. ...... ADD  IN  CONTRIBUTIONS  DUE  TO  SINGULARITIES 

I  20.  NJ«244(N— I-l) 

i30.  OO  220  K  *  1*  NJ 

>31.  JVK  =  JV(I  4  I.K)  /  2 

\  32.  Jmk  =  Jw (  I  4  i.K)  /  2 

;  33.  SVK  =  SVC  I  4  I .K) 

J34.  S*K  =  SM<  I  4  I  *K) 

135.  VCI.l  -  JVK. 2)  =  VCI.l  -  J  VK • 2 )  4  AM  4  SVK  *  XM(  I  •  1 i 

t  36.  MCI. I  -  JVK.2J  *  MCI.l  -  JVK. 2)  4  AP  4  SVK  *  XP<  I  .  1  J 

137.  NS  s  NTGP  -14  JVK 

)3B.  NF  s  2  -  JVK 

139.  DO  160  K I  =s  NF.  NTGP 

lAO.  DO  150  M  *  I.  2 

(41.  V(I.KI.M)  a  V(l.KI.M)  ♦  AM  *  SVK  4  XMCI.KI  4  JVK) 

(42.  M<  1 . K I .M )  -  *<I.K1.M>  4  AP  *  SVK  *  XP(  1  .KI  4  JVK) 

(43.  150  CONTINUE 

144.  160  CONTINUE 

145.  IF  (NS  -  Nl)  180.  17  0.  170 

146.  170  KKI  =  NP1  -  JVK 

147.  V(I.NP1  -  JVK. 2)  =  (V(I.NPl  -  JVK. 2)  -  AM  •  SVK  *  XM(l.NPl)) 

148.  1  4  AM  *  SVK  *  XMCCI) 

149.  MCI. NP 1  -  JVK. 2)  *  (MC1.NP1  -  JVK. 2)  —  AP  4  SVK  4  XPCI.NPIJ) 

150.  1  4  AP  4  SVK  *  XPCCl) 

151.  180  NS  «  NTGP  -  1  4  X( I  4  1)  -  JMK 

152.  NF  *  JNK  -  X(  1  4  1  )  4  2 

153.  VCI.NF  -  1.2)  =  VCI.NF  -  1*2)  ♦  AM  4  SMK  *  YMCI.l) 

154.  «<I. NF  -  1.2)  •  MCI. NF  -  1.2)  4  AP  4  SMK  4  YP(l.l) 

155.  DO  200  K 1  =  NF.  NTGP 

l5o.  KKI  =  KI  4  XCI  4  1)  -  JMK 

157.  00  190  M  -  1.  2 

158.  VCI.KI.M)  =  VCI.KI.M)  4  AM  4  SVK  4  YM(I.KKI) 

159.  VCI.KI.M)  =  4C1.KI.M)  4  AP  4  SMK  4  YP(I.KKI) 

160.  190  CONTINUE 

161.  200  CONTINUE 

162.  IF  (NS  -  Nl)  220.  210.  210 

163.  210  KKI  *  JVK  -  X(I)  4  1 

164.  V ( I • KK 1.2)  *(VCI.KKI.2)  -  AM4SMK4VMC I . NP l ) )4  AM4 SmK4 YMC ( I > 

165.  MCI. KKI. 2)  *<M(I.KKI.2)  -  AP4 S WK4 YP ( I . NP 1 ) )  4 AP4S »K*Y PC ( I i 

166.  220  CONTINUE 

167.  230  CONTINUE 

163.  240  CONTINUE 

169.  C....... CALCULATE  SCATTERING  KERNELS  MS  ANO  VS 

170.  CP  *  (C(l)  4  l.)  /  2. 

171.  CM  *  <C( 1)  -  1.)  /  2. 

172.  DO  260  K  =  1 .  NTGP 

173.  DO  250  M  *  1.  2 

174.  VS(K.M)  *  -  CM  4  V(l.K.M)  4  CP  4  M(l.K.M) 

175.  WSCK.M)  =  CP  *  V(l.K.M)  -  CM  4  M(l.K.M) 

176.  250  CONTINUE 

>77.  260  CONTINUE 

>78.  CALL  PR  NT ( 2 ) 

179.  RETURN 

130.  =ND 
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SUBROUTINE  FIELD 

•  ••••THIS  SUBROUTINE  COMPUTES  SCATTERED  ANO  INTERNAL  FIELDS 


c . -..THROUGH  use  of  THE  SCATTERING  KERNELS  OBTAINED  IN 

- . . ...SUBROUTINE  DELTA  AND  THE  R I £MANN  FUNCTI ONS  OBTAINED 

C . ••••IN  5U6RQUT  I  Nc  KERNEL 


INTEGER  X 

DIMENSION  UT12001  )  •  **(  129 i  • S ( 3  2)  * J  J( 32 ) 

COMMON  / BLOCK 1 /H*  NL • NTGP • NPTLt 5) • I FRO# ITM *  IT » I SCA T • I OPT 
COMMON  /BL0CK2/A( 33) #6(33) «AP<  33) *B°<33) • DTC33 1 • AHC 33 )  *BH(33). 

1  APH ( 33 )  . BPH ( 33 ) •DTHC33 )  . AI  ( 331 #81 ( 33) tOT| (33) . 

2  OP(33) «E8 I (33). AR( 5 ).AL (5) .8R (5). 8l( 5) 

COMMON  /9L0CK3/X ( 6 ) »C ( 6). DPI  133) .OPBIC33) .EAMBC5.33) .EAPBC5.33) . 

1  QOLDC 33 ) • ONEm (  33  )  . XP( 5# 129 ) . XM C 5. 1 29 > . YP< 5. 129)  • 

2  YMC  5.  129)  •  XPCl  5  )  *  XMC(  5)  *  YPC(  5)  •  YMCC  5)  .APBUS)  * 

3  AMB 1(5) #P( 5.2. 33.33) 

COMMON  /BLOCK4/JV(  5.16) .UMCS.16) .SVC  5. 16 ) .SMC 3# I  6 ) • VC  5 #200 1 #2) • 

1  ■( 5. 129.2) •  VS( 1 29.2) .*5(129.2) 

COMMON  /BLOCKS/ BUFF ( 2201 1 • UI ( 2001 ) .UC 200 1 i .UR I  200 1 > 

EQUIVALENCE  ( BUFF ( 130 ) .UT C 1 ) ) 

REALMS  SI .S2.S3.S4.S5.S6 
DO  10  I  *  1.  2200 
BUFF ( I )  *  0. 

10  CONTINUE 

. •••COMPUTE  TRANSMITTED  FI EU) 

km  a  2  MM  (NL  -  1  ) 

CP  *  (Cl  1 >  ♦  1.)  /  2. 

CM  »  ( C( 1 )  -  1.1  /  2. 

OO  20  K  *  1  t  KM 
JJCK)  *  JV(l.K)  /  2. 

JJCKM  ♦  K)  a  -  JttCl.K)  /  2. 

5(  K)  *  SV(l.K)  *  CP 
S(KM  ♦  K)  m  -  Stf(l.K)  M  CM 
20  CONTINUE 

00  30  M  «  2.  NTGP 
WW(M)  *  CMS(M.l)  ♦  MSCM.21)  *  H 
30  CONTINUE 

WMC  *2.  •  H  •  MS ( NTGP • 2 ) 

D 1 V  *  SCI)  ♦MM  VS( 1 .2) 

UT( 1 )  x  0. 

UT(2)  *  UI ( 2)  /  DIV 


4  l  .> 

42. 

43. 
44  . 
45- 

46. 

47. 
43. 

49. 

50. 
51  . 

52. 

53. 

54. 

55. 

5  6  . 
57  . 
5  3* 
59. 
60  . 


IM  *  2  ♦♦  NL 
00  60  K  *  2.  NTGP 
KP  *  K  ♦  1 

51  *  0.00 

52  a  0 .DO 

DO  40  I  «  2 •  K 

51  *  SI  ♦  MM(I)  4  UT(2  ♦  K  -  I) 

40  CONTINUE 

00  SO  I  a  2.  IM 

52  «  S2  ♦  SCI)  M  UTCKP  ♦  JJCI)) 

50  CONTINUE 

UT(KP)  *  (UI(KP)  -  SI  -  S2I  /  OlV 
60  CONTINUE 

KF  a  NTGP  ♦  1 

52  a  0.00 

ITMI  a  ITM  ♦  1 

DO  90  K  a  KF.  ITMI 

KP  a  K  ♦  l 

SI  a  0 .DO 

DO  70  I  *  2.  NTGP 
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ol  • 

si  =  SI  ♦  «  uT  (  2  «■  K.  -  I) 

62  • 

70 

CONT INUc 

63* 

s2  =  52  ♦  UTCKP  -  NTGP) 

S3  =  0.00 

65* 

DO  30  I  =  2.  IM 

66, 

S3  =  S3  +  S ( I )  *  UT( KP  ♦ 

67. 

60 

CC NT INUE 

oB. 

UT(KP)  =  (UK  KP)  -  S 1  -  C  *  S2 

69. 

90 

CONTINUE 

70. 

C  .  .  .  • 

71  . 

DO  91  K  s  1,  KM 

72. 

S(K)  =  -  CM  *  S(K)  /  CP 

73. 

S(KM  ♦  K)  ss  -  CP  ♦  S  (  KM  ♦  K)  /  CM 

74. 

91 

CONTINUE 

75. 

DO  92  M  =  2,  NT OP 

76. 

■"»<*)  =  (VSCM.l)  «■  VS ( M » 2 )  )  *  H 

77. 

92 

CONTI NUE 

73. 

**C  3  2.  *  H  *  VS1NTGP.2) 

79. 

S(l>  =  5(1)  4  H  *  V5< 1.21 

80. 

UR  (  1 )  =  0. 

01  • 

UR  12)  =  S(i)  *  UT (2) 

e2. 

DO  95  K  =  2.  NTGP 

S3. 

KP  =  K  ♦  i 

e4. 

SI  *  0.00 

es. 

S2  =  0.00 

86. 

OO  93  I  =  2*  K 

87. 

SI  *  SI  ♦  MW< I )  *  UT ( 2  ♦  K  —  I ) 

98. 

93 

CONTINUE 

89. 

DO  94  I  —  1«  IM 

90. 

S2  =  S2  +  S(I)  *  UTCKP  JJ(IJ) 

91  . 

94 

CONTINUE 

92. 

URtKP)  =  51  ♦  S2 

93. 

95 

CONTINUE 

94. 

S2  *  0.00 

95. 

DO  98  K  =  KF.  ITM 

96. 

KP  »  K  ♦  1 

97. 

SI  =  0.00 

98. 

DC  96  1*2.  NTGP 

99. 

SI  *  SI  ♦  fcw(I)  *  UT ( 2  ♦  K  -  I) 

1  00. 

96 

CONTINUE 

101. 

S2  *  S2  ♦  UTCKP  -  NTGP) 

102. 

S3  =  0.00 

1  03. 

DO  97  I  *  l.  IM 

1  04. 

33  =  S3  4  SCI)  *  UTCKP  ♦  J JC I ) ) 

1  05. 

97 

CONTINUE 

1  06. 

UR  C  KP )  *  SI  +  52  *  WWC  ♦  S3 

1  07. 

98 

CONTINUE 

1  03. 

DO  100  K  =  1.  ITMl 

1  09. 

V(NL.K.l)  *  UTCK) 

1  10  . 

V C  NL • K .2 )  3  0. 

111. 

100 

CONTINUE 

1  12. 

CALO  PRNTC3) 

i  13. 

C . 

1  14. 

IF  (NL.EQ.l)  GO  TO  160 

115. 

NT  D  *  NL  —  1 

1  16. 

OO  150  IL  *  1 ,  NTD 

:  17. 

NCL  *  NC  -  IU  ♦  1 

1  13. 

N»L  *  NCL  ♦  l 

119. 

NNL  *  NCU  -  1 

120. 

CV1  *  ^BI(NCL)  *  <C(N  PL)  ♦  1.) 
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L  2  1  • 
122. 
123* 
12*. 

I  25 • 

1  26. 
127. 
126. 

129. 

130. 
131  • 

132. 

133. 

134. 

135. 

136. 

137. 

138. 

139. 
1  40  • 
l  41  . 

142. 

143. 
I  44  . 
145. 
1  46* 

147. 

148. 

149. 

150. 
IS!  • 
152. 
I  S3. 

154. 

155. 

1 56. 
157. 
1  56* 

159. 

160. 
161  . 
162. 
1  63. 
I  64  • 

165. 

166. 
167. 
t  68. 

169. 

170. 

171. 

172. 

173. 

174. 
I  75. 
I  76. 
1  77. 
I  79. 

I7?. 
I  90. 


CKl  *  -  APBKNai  •  (C(NPL)  -  !•)  /  2. 

CV2  *  -  A  MB  I  (NCi.)  *  ( C  (  NPL )  -  1  .  )  /  2. 

C»  2  —  APB  I  (NCI.)  *  (C(nPl)  ♦  U)  /  2* 

CV3  =  h  *  AMBKNCL)  /  2, 

C*3  *  h  *  APBI(NCL)  /  2. 

CVA  =  H  *  XMC(NCL)  •  PMBKNCl)  /  A. 

CVA  *  W  *  XPCt  NCL>  •  APB  I  ( NCI.)  /  A. 

CV  5  »  h  *  VMCCNCL)  *  AMBMNCU  /  A. 

CVS  *  H  *  YPC  ( NCL )  •  APBKNCU  /  A. 

NPl  *  NPTL(NCL) 

NPL  *  NPL  -  1 
V(  NNL.  M)  *  0* 

V ( NNL • l .2)  *  0. 

00  120  <  »  2*  NPl 

51  *  -  XM(NCL.l)  •  V(NPL.K.l)  ✓  2. 

52  »  -  YM(NCL.l)  •  V(NPL.K.2>  /  2. 

53  *  —  XP( NCL  > 1 )  A  V(NPL.K.l)  /  2. 

SA  X  -  YP(NCL.l)  A  VC  NPL.K  >2 )  /  2* 

KM1  *  K  -  1 

00  110  J  »  1.  KM1 
KP  A  K  A  l  —  J 

51  *  SI  A  XM(NCL.J)  A  V ( NPL * KP *  1 ) 

52  *  52  »  YM(NCL.J)  a  V( NPL.KP, 2) 

53  *  S3  A  XP(NCL.J)  A  V(NPL.KP.l) 

SA  *  SA  ♦  YP(NCL.J)  *  V(NPL.KP*2) 

110  CONTINUE 

V ( NNL*  K. 1 )  *  CV1  A  V(NPL.K.I)  ♦  CV2  •  V(NPL*K.2)  A  CV3  • 

1  (SI  A  S2 ) 

V( NNL*  K»  2)  •  C«3  *  (S3  A  SA) 

120  CONTINUE 

55  »  0.00 

56  »  0.00 

NP2  »  NPl  A  1 

NM 1  *  NPl  -  1 

OO  1 AO  K  «  NP2.  I  TNI 
KM  «  K  -  NM1 

51  «  (XM(NCL.l)  •  V(NPL.K.I)  A  XM(NCL.NPl)  A  V(NPL.KM.D)  / 

52  *  (YM(NCL.l)  *  V( NPL.K. 2)  A  YM(NCL.NPt)  •  V(NPL.KM.2I)  / 

53  s  (XP(NCL.l)  A  V( NPL.K. I)  A  XP(NCL.NPl)  A  V(NPL.KM.l))  / 

SA  a  (YP(NCL.l)  *  V( NPL.K*  2)  A  YP(NCL.NPl)  •  V(NPL.KM.2))  / 

DO  130  J  a  2,  NM 1 

KP  a  K  A  1  -  J 

51  *  SI  A  XM(NCL.J)  *  V(NPL.KP.l) 

52  •  S2  A  YM(NCL.J)  A  V ( NPL.KP .2 ) 

53  «  S3  A  XP(NCL.J)  •  V ( NPL.KP . 1 ) 

SA  a  SA  A  YP(NCL.J)  «  V( NPL.KP .2 ) 

130  CONTINUE 

55  «  S5  A  V(NPL.KM.l)  A  V(NPL.K-NPl.t) 

56  «  S6  A  V( NPL .KM .2)  A  V(NPL.K-MP1.2) 

V(NNL.K.l)  a  CV1  •  V(NPL.X.l)  A  CV2  •  V( NPL.K. 2)  A 

1  CV3  «  (SI  A  S2)  A  CVA  *  SS  A  CV5  •  S6 

V( NNL.K.2)  ■  C«1  *  V(NPL.KM.l)  A . C«2  •  VCNPL.KM.2)  A 

1  C«3  •  (S3  A  SA)  A  CM  *  SS  a  CVS  A  $6 

1  AO  CONTINUE 
ISO  CONTINUE 
160  CONTINUE 

M2  «  M  *  2. 

JC  170  I  a  l,  ITM1 
UT( I >  •  1  * 

170  CONTINUE 
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M  M  h  h 


1  31  . 

152. 
J  S3  • 

I  04. 

I  S3, 
i  96 * 
I  57. 
;  as. 

IS9. 
190. 
I  91  . 

192. 

193. 
I  94. 

195. 

196. 

197. 
193. 

1  99. 
200. 
201  • 
202. 

203. 

204. 

205. 

206. 
207. 
203. 

209. 

210. 
211. 
212. 

213. 

214. 

215. 

216. 

217. 

218. 

219. 

220. 
221  • 
222. 

223. 

224. 

225. 

2  2e>. 
227. 

226. 

229. 

230. 

231. 

232. 

233. 

234. 

235. 
23o. 
237. 
233. 
229. 
240. 


IF  (I0PT.5Q.1)  SO  TCJ  300 

C. ........ .CALCULATE  FIELD  INSIDE  MEDIUM  AT  TIMES  T  a  IT  *  H  *  2 

DO  220  IT  »  IFRQ.  I T M  *  IFRQ 
IN  =  NTGP 

IJ=  IT  4  1—  < NTGP  —  l I  /  2 

U (  IN)  =  V ( NL • I J • 1  J  *  UT(IJ) 

DO  210  K1  =  1 .  NL 

NCL  =  NL  -  Kl  4  1 
NPL  *  NCL  4  1 
NP1  *  NPTL(NCL) 

NP2  s  NP1  M 
NM  l  a  NP  i  -  1 

CM  =  C(NPL)  -  l. 

CP  *  C(NPL)  4  1. 

NIT  s  <NPl  -  U  /  2 
00  200  JC  *  1.  NIT 
IN  *  IN  -  2 

M  *  NIT  -  JC 

M2  *  M  4  2 

IJ  *  IT  -  M  4  1  -  X(NCL)  /  2 

IJM  <=  I J  —  NM  1  4  M2 

51  *  ( V ( NCL  t I J • 1  )  *  PINCL.1. M2  4  1.1)  *  UT (  I J J  ♦ 

1  V(NCL.IJM.I)  4  P(NCL. 1.NP1.NP1  -  M2)  *  JT(IJM))  /  2. 

52  *  ( VtNCL.I JM.2)  4  P < NCL .  2  •  M  2  4  1.1)  4  UT(IJM)  4 

1  VINCL.IJ.2)  4  P(NCL.2.NP1 .NP1  -  M2)  4  UT(IJ))  /  2. 

IF  ( I J.LE.l )  GO  TO  190 
JM  *  2  4  JC 

IF  ( I JM.LE. 0 i  JM  *  I J  —  l 
DO  180  J  *  2*  JM 
M2J  *  M2  ♦  J 

51  a  Si  4  V ( NCL .  I  J  41-  J.l)  4  P( NCL. I .M2J. J) 

52  *  S2  4  V ( NCL . I J  41-  J.2)  4  P( NCL. 2.NP2  -  J. NP2  -  M2 J J 

180  CONTINUE 

190  CONTINUE 

MI  2  *  M2  4  1 

U(  I N )  a  EAMBi  NCL.MI2)  4  (CP  4  V(NCL.U.l)  4  UT(IJ)  - 
1  CM  4  V ( NCL . 1 J . 2)  4  UT(IJ)  4  H2  4  SI) 

U( IN)  a  U( IN)  4  EAPBCNCLt MI2)  4  (CP  4  V ( NCL. I JM. 2 )  4  UT(IJM) 
1  CM  4  V(NCL.IJM.I)  4  UT ( IJM)  4  H2  4  S2 ) 

200  CONTINUE 
210  CONTINUE 

CALX  PRNT(4> 

220  CONTINUE 
RETURN 

300  CONTINUE 

. . CALCULATE  FIELD  INSIDE  MEDIUM  AT  POSITION  IX 

RE  AO  (  5.1000)  NOB  PT  5 

00  390  KK  =  1  •  NOBPTS 
READ  (5.1000)  IX 

IF  (IX  -  NTGP  4  l)  330.  310.  330 
310  CONTINUE 

JM AX  a  1  4  (ITM  -  l  -  (IX  /  2)  )  /  IFRQ 

DO  320  J  a  1,  JM A  X 

U(J)  a  VCNL.IFRQ  4  (J  -  l)  4  1.  1) 

220  CONTINUE 
GO  TC  380 
230  CONTINUE 

NLP  a  NL  4  l 

00  340  1=2*  NLP 

IF  (IX.LE.XCl))  SO  TO  350 
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-41  . 
242  • 
243. 
?  44  • 
245. 
246  • 
247. 
?45. 

249. 

250. 
2Si  • 

252. 

253. 
2  54. 

255. 

256. 
2  57, 
253. 

259. 

260. 
*61  • 
262. 

263. 

264. 

265. 

266. 

267. 

268. 

269. 
2  70. 
271  . 
272. 
’73. 

274. 

275. 

276. 

277. 

270. 
2  79. 
280. 
231. 
262. 
283. 
2  34. 
’35. 
2  86. 


340  CONTINUE 

*«ITE  (6.2000)  IX 
RETURN 

350  CONTINUE 

NCL  *  I  -  1 

NPL  *  NCL.  ♦  l 
NPi  *  NPTL(NCL) 

NP2  *  NP1  ♦  1 
NM|  *  NPI  -  1 
CM  *  C(NPL)  -  l. 

CP  *  Cl NPL)  ♦  1. 

NIT  «  (NPI  -11/2 
M  «  (IX  -  X(NCL) )  /  2 
M2  s  M  •  2 
MI  2  *  M2  ♦  1 
JC  *  NIT  -  m 
U(  1  )  *  0  • 

IMAX  «  1  ♦  (ITM  -  (IX  /  2)  )  /  I  PRO 

00  370  I  *  2.  IMAX 
I J  »  IFRO  ♦  ( I  -  1 )  ♦  1 
IJM  *  1J  -  NMi  ♦  M2 

51  *  ( V( NCL • I J  t  X )  •  P ( NCL. 1. M2  +  1*1)  *  UT(1J)  ♦ 

1  V(NCL.IJM.i)  4  P( NCL. 1 .NPI. NPI  -  M2)  4  UTCIJM))  /  2. 

52  *  ( V( NCL. I JM.2)  4  PCNCL.2.M2  ♦  1.1)  6  UT(lJM)  4 

1  VINCL.IJ.2)  4  P(rCL«2.NPl .NPI  -  M2)  4  UT(IJ))  /  2. 

JM  *  2.  4  JC 

IF  (IJM.LE.O)  JM  »  IJ  -  1 
00  360  J  *  2«  JM 
M2  J  9  M2  ♦  J 

51  *  SI  ♦  V1NCL.IJ  ♦  I  —  J<>  )  4  PC  NCL*  1.M2J.J) 

52  *  S2  ♦  V(NCL.XJ  ♦  1  -  J o2 )  *  PCNCL.2.NP2  -  J. NP2  -  M2Ji 
360  CONTINUE 

MI  2  *  M2  ♦  1 

UCI)  «  EAMBC  NCL.  Ml  2)  4  C  CP  4  V(NCL.lJ.I)  4  UT(IJ)  - 
I  CM  4  V ( NCL* I J • 2 )  «  Ut ( I J )  ♦  H2  •  SI) 

U(I)  «  UCI)  4  EAPBCNCL.MI2)  4  (CP  4  V ( NCL • I JM. 2)  4  UT(IJM) 

l  CM  4  V  C  NCL • I JM » l )  4  UT ( I JM)  ♦  H2  4  S2 ) 

370  CONTINUE 
380  CONTINUE 
IT  *  IX 
CALL  PRNT(S) 

390  CONTINUE 
RETURN 

1000  FORMAT  CI5) 

2000  FORMAT  ( /// .IX. 9HHEY .  IX  «*I5.1AH  IS  TOO  LAROE) 

END 
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VIII.  SETTING  DP  SOME  PROGRAMS 


Consider  a  three-layer  medium  of  depth  3  cm,  with  each  layer  beinq  I  cm 
deep.  Assume  that  the  permittivity  and  conductivity  are  constant  on  each 
layer  and  are  given  by  the  values  shown  in  Figure  B-3. 

N  is  chosen  to  be  128.  This  is  not  for  purposes  of  accuracy  (N=64  yields 
identical  results)  but  rather  for  obtaining  more  grid  points,  resulting  in 
smoother  plots. 

Turning  to  the  computation  of  the  X(i)'s,  note  first  that 
t=(0.01)(/30  +  /135  +  /60)/c 

where  c  is  the  speed  of  light  in  free  space.  Then  X(1)=0  and 
x1=/30/(/310=/T35  +  /60)=X(2)/128. 

Thus,  X(2)=28.22  •••  and  so,  upon  rounding 
X(2)=28. 

This  is  equivalent  to  using 

e1(z)=29.4eQ  rather  then  et(z)=30c0. 

Now 

x2=(/TO‘  +  /l35)/(/29.4  +  /T3S-  +  /S5) 

=X(3)/128 

and  so 


X(3)=88 
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files  for  a  three-layer  medium. 


and 


X(4)=l?8. 


Since 


X(4)-X(3)>32 

X(3)-X(2)>32 

two  additional  layers  need  to  be  inserted  into  the  medium.  (Note  that 
N=64,  no  additional  layers  are  needed.)  Hence,  set 

X(1)=0 

X(2)=28 

X(3)=60 

X(4)=88 

X  ( 5 ) = 1 08 

X(6)=128. 


Now 


C(l)=/2974 

C(2)=/135/29.4 

C(3)=1.0 

C(4)=/WM 

C(5)=l  .0 

C(6)=/T7M. 

The  only  nonzero  functions  to  be  entered  into  subroutine  EVAL  are 

Bl(x)=-(1.0)t/(29.4e0) 

B2(x)=B3(x)=-(O.l)*/(135e0) 

B4(x)=B5(x)=-(0.5)t/(60e0) 
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A  1000-MHz  sinusoidal  Incident  field  of  duration  20  ns  is  used.  Thus, 

E1  (-t)=s1n(u>t)»H  (20»10_9-t) 

C 

9 

where  u>=2ifl0  and  Hg  Is  the  Heaviside  function.  Since  T=t/t,  It  follows 
that  the  incident  field  to  be  used  In  INCWAV  is 

UINC(T)=s1n(-u>lT) . 

Now  the  t  step  size  is 

Ar=2H=l/64 

so  it  follows  that  the  t  step  size  is 
At=*/64=. 01291  ns. 

The  period  of  the  incident  signal  is  1  ns  which  corresponds  to  approximately 
77.5  At.  In  light  of  this  the  incident  signal  was  modified  so  that  its  period 
was  an  even  multiple  of  At.  (This  modification  Is  not  at  all  necessary  and 
was  performed  solely  to  make  numbers  turn  out  "nice.")  The  period  of  E1  was 
lengthened  to  78  At  which  corresponds  to  a  993-MHz  signal.  Then  UINC  can  be 
written 


UINC(x)*s1n(-it«T/(78»H) ) . 

(See  Section  IX  for  a  source  listlnq  of  INCWAV  for  this  particular  field.) 

The  duration  of  this  new  Incident  field  was  lengthened  to  20.13951  ns  so 
that  20  cycles  would  still  be  considered.  Thus,  ITM  should  be  at  least  1560 
so  that  signal  cutoff  Is  observed.  It  was  decided  to  print  out  the  entire 
interior  field  every  half  period,  so  IFRQ  was  set  equal  to  39  and  ITM  was  cho¬ 
sen  to  be  1950.  Setting  ISCAT*1  and  I0PT*0,  the  scattered  and  Internal  fields 
were  put  onto  disk  storage.  These  fields  were  then  graphed  using  the  routine 
given  in  Section  X. 
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The  above  results  indicate  that  steady  state  was  essentially  attained 

after  six  periods  of  the  incident  field;  i.e.,  after  468  At  intervals.  It  was 

then  decided  to  search  for  the  maximum  intensity  occurring  at  each  grid  point 
during  the  passing  of  the  transient  field  (first  6  cycles)  and  during  steady 
state.  mus,  HATS  was  run  a  second  time,  with  IFRQ=1,  ITM=468,  ISCAT=0,  and 
I0PT=0  and  with  the  modified  subroutine  PRNT  given  in  Section  XI.  To  deter¬ 
mine  the  maximum  steady  state  Intensities,  HATS  was  run  with  IFRQ=780  and 

ITM=858.  Upon  the  first  call  to  PRNT  usinq  PRNT(4),  the  value  of  I  FRO  was 

changed  to  1  so  that  the  entire  period  was  examined  in  as  much  detail  as  pos¬ 
sible.  Results  from  these  runs  were  plotted  on  the  same  set  of  axes. 

From  these  plots  it  was  observed  that  two  points  in  the  medium  (one  in 
the  first  layer,  one  in  the  third)  had  a  considerably  stronger  transient  field 
than  steady  field.  Thus,  a  fourth  run  of  HATS  was  made  lookinq  at  the 
detailed  time  behavior  at  these  points  as  well  as  at  a  point  in  the  second 
layer  which  seemed  to  have  almost  no  difference  in  transient  vs.  steady  state 
intensity.  Thus,  three  more  plots  were  generated  with  the  parameters  used  in 
HATS  being  IFRQ=3;  ITM=1950;  ISCAT=0;  I0PT=1;  N0BPTS=3;  and  IX=20,  32,  92. 

For  the  sake  of  interest  this  entire  procedure  was  pursued  a  second  time 
for  a  nonconducting  medium.  This  was  achieved  by  settinq  Bl(x)=B2(z)=... 
=B5(x)=0. 

Finally,  so  that  the  internal  reflections  in  the  medium  in  the  non¬ 
conducting  case  could  be  observed,  a  triangular  incident  spike  field  was  used: 

IJI  (1  )=0 
UI (2)=0.5 
UI(3)=1.0 
UI  (4)=0.5 
UI(5)=UI(6)=***=0. 

Plots  of  all  output  appear  in  Section  XII. 
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IX.  SOURCE  LISTING  OF  INCWAV  FOR  FIELD  OF  FINITE  DURATION 


t.  5U8HCUT I  imc  JNCWAV 

2.  C.  .........  I NSERT  FUNCTION  STATEMENT  FOB  INCIDENT  FIELD 

3.  JINC(X)  *  S I N (  -PI  »  X  /  (7B.  *  N>  * 

4*  C. ........ . 

5.  CCMMCiN  /BLOCK.  l/H.NL.NTuP.NPTL(  5)  .  IF«0.  ITM,  IT  .  I  SC  AT  •  IOPT 

fc.  COMMON  /3LUCKS/BUFF( 2201 ) . UI ( 2001 > ,J(200 l I •  UH  ( 200 l > 

7.  PI  »  3.1*15927 

8.  H2*M*2. 

9.  ITM1  *  ITM  ♦  I 

10.  DC  10  J  *  l.  ITM1 

11.  JI C J)  »  0. 

12.  10  CONTINUE 

13.  DO  20  J  *  2.  39 

14.  X  *  -  (J  -  1)  *  H2 

15.  UI  (  J)  *  UINCCX) 

16.  UI ( 39  ♦  J)  *  -  UII JJ 

17.  20  CONTINUE 

13.  DC  40  J  *  2.  20 

19.  JC  *  I J  -  1  )  *  78 

20.  DO  30  I  *  2*  78 

21.  UI C JC  ♦  I >  »  UIC I  ) 

22.  30  CONTINUE 

23.  4 0  CONTINLi 

2*.  RETURN 

25.  eNO 


X.  SOURCE  LISTING  OF  A  PLOTTING  ROUTINE 


1 . 

2. 
3. 
A. 


24. 

25. 

26. 


DIMENSION 
31  MENS  ION 


READ 


X( 193). Al (21 00 1.A2 165). A3( 22001.0(193) 

X0(2).X1(2).X2(2).X3(2).A4(2).U0(2).U1(2).J2(2). 

03(2 ) .U4(2) 

(8.2000)  N.  NL.  IFRQ.  ITM 


5. 

N  *  N  /  2 

6  • 

H2  —  1  •  /  N 

7. 

NP1  *  N  ♦  1 

9. 

N2  *  N  *  2 

9m 

N2  1  *  N2  ♦ 

1 

10. 

MN  *  N  *  3 

♦  1 

11* 

00  1 0  I  * 

Is  NN 

1  2. 

X(l)  *  *1  . 

♦  <  I  - 

13. 

10 

CONTINUE 

14. 

DO  20  1  * 

1  t  N 

15* 

Al  (1  )  *  0* 

16* 

20 

CONTINUE 

17. 

DO  30  I  * 

1  •  N2 

16* 

A3( 1 )  »  0. 

19. 

30 

CONTINUE 

20. 

ITM1  »  ITM 

♦  1 

11. 

L  •  ITMi  ♦ 

N2 

22. 

RE AO  (8.1500)  (A 

22. 

L  *  ITNI  4 

N 

*  H2 


1 


N21.  L) 


READ  (8.1500) 
C.« ........ EST AflLISH 

. . LAVERS 


(Aid).  I  •  NP1  .  L) 
ARRAYS  TO  BE  USEO  IN 
ANO  HORIZONTAL  AXIS 


GRAPH! NO  BOUNDARIES  OP 
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II 


27 .  XOil)  a  0. 

28.  X0C2)  *  X0<  1  ) 

29.  XI (i)  a  .21875 

30.  XI <  2 )  a  Xl<  1) 

31.  X2 ( 1 )  *  .6875 

32.  X2(2J  a  X2( 1 ) 

33.  X3 ( 1 )  a  x.O 

34.  X3 C  2 )  a  X3(l) 

35.  X4C1 )  a  -1. 

36.  X*<2)  a  2 . 

37.  U0< 1 )  =-1.2 

38.  UO (21  a  1.2 

39.  Ul  <  1  )  =  -.75 

40.  U 1 ( 2 )  a  *75 

4  1.  U4  (  1  )  a  0 . 

42.  U4<2)  a  0. 

43.  DO  100  IT  a  I FR  3 «  ITM  *  IFRQ 

44.  READ  ( 6*6000)  K 

45.  READ  ( 8 1  1500)  (42(1).  1  =  1*  NPl ) 

46.  OC  40  I  =  1.  NPl 

47.  U<  I  )  =  A1  (  IT  ♦  I  ) 

48.  U(N2  «■  I)  a  A3  (  NP  l  ♦  IT  -  I  ♦  1) 

49.  40  CONTINUE 

50.  DO  50  I  =  l,  NPl 

51  .  U<  N  4  I  )  =  A  2  (  I) 

52.  50  CONTINUE 

53.  C... . THE  ARRAY  U  NO*  CONTAINS  THE  REFLECTED.  INTERNAL  AN 

54.  C. ............ .TRANSMI TTEO  FIELDS  AT  TIME  TAU  =  IT  *  2  *  H 

55.  C.  . ..PLOT  ARRAY  J 

56.  CAlL  GRAPH  ( MN «  X.  U»  0.  2  *  6.01.  4. 01.  . 50t— l.. .60.  —  1.2 

57.  1  ’NORMALIZED  DEPTH:*,  ’NORMALIZED  INTENSITY’, 

5a.  2 

59.  . . PLOT  VERTICAL  LINES  SHOEING  POSITION  OF  LAYERS 

60.  CALL  GRAPHS  (2.  XOi  UO,  0.  4,  •;•) 

61.  CALL  3AAJHS  (2.  XI,  Ul,  0.  4.  •;•) 

62.  Cm  Li-  GR  APHa  (2,  X  2  .  Ul.  0,  4.  *»*) 

63.  CALL  GRAPHS  (2,  X3,  UO,  0.  4,  •;*) 

64.  C-... . »LOT  HORIZONTAL  AXIS  AT  ZERO  INTENSITY 

65.  CAll  GRAPHS  (2.  X4,  OU ,  0,  4,  *;•) 

66.  100  CONTINUE 

67.  STOP 

69.  1500  FORMAT  (8E10.3) 

69.  2000  FORMAT  (415) 

70.  6000  FORMAT  (110) 

71.  END 


XI.  SEARCHING  FOR  MAXIMUM  INTENSITY 


The  following  source  listing  is  a  modification  of  subroutine  PRNT  used  to  search 
the  program  output  for  maximum  intensity. 


1.  SUBROUTINE  PRNT(K) 

2.  INTEGER  X 

3.  COMMON  /8l0C<I,/’H.NL,NTGP.NPTL(  5)  .IFRQ.ITM.IT.  ISCAT.iOPT 

A*  COMMON  /BL0CK2/A{  33>.a(33).AP(33).BP(33).DT(33).AH(33).BH(33). 

5.  1  A=>m(33J  .3PH(33J  .0TH(33)  .AI  (331  .81  (33)  .OTI  (33)  , 

6.  2  □»{ 33 ) . E3I ( 33 ) . AR(5) . AL( S) *8R(S) . 8L ( 5 ) 

7.  COMMON  /  3LCCK3/  X  (  6  )  .  C  (  6  )  .  OP  I  (  33  )  .OPS  I  (33  )  .EAMB<5 ,33)  . EAP6(  5.33  )  . 

3.  1  OOLD( 23 ) ,QNEW (33) . XP( 5. 129 ) ,XM( 5, 129) , YP( 5. 129) . 

9.  2  YM(  5,  129)  .  XPC(  5  )  •  XMC(  S)  .  YPC<  51  .  YMC(  5)  •  APBK  5)  . 

10.  3  AMBI( 5) ,P< 5.2,33.33) 

11  .  COMMON  /BL0CK4/JV15.  16) , JW(5,16> .SV(5. 16) .SW( 5.16) .VC  5.2001.2)  . 

12.  1  w<5.129,2).VS(129.2).WS(129.2) 

13.  COMMON  /BLOCK5/8UFF(2201) ,UI (2001 ) .U( 2001 ).UR( 2001 ) 

1*.  DIMENSION  £MAX(129).ITI ME( 129) 

15.  IF  IK  -  2)  1.  2,  50 

16.  1  CONTINUE 

17.  C.. ..... ...PRINT  OUT  INPUT  DATA  FROM  LINES  1.  2.  3 

13.  N  =  NT GP  -1 

19.  WRITE  (6,2000)  N,  NL.  IFRO.  ITM .  ISCAT.  IOPT.  (X(I).  I  =  1,  NL) 

20.  WRITE  (6.3000)  (C(I).I  *  1.  NL) 

21.  C  WRITE  (3.2000)  N.  NL .  I FRO ,  ITM 

22.  RETURN 

23.  2  CONTINUE 

2*.  RETURN 

25.  50  CONTINUE 

26.  IF  (K  -  4)  3.  4,  60 

27.  3  CONTINUE 

25.  C. ...... ...PRINT  INCIDENT.  REFLECTED  AND  TRANSMITTED  FIELDS  IF  1  SCAT  = 

29.  IF  (ISCAT  .EO.  0)  RETURN 

30.  ITM1  =  ITM  ♦  1 

31.  WRITE  (6.5000) 

32.  WRITE  (6.1000)  (UI(I).I  =  1,  ITM1 ) 

33.  WRITE  (6.5100) 

34.  WRITE  (6.1000)  (UR(l).I  *  1,  ITM1 ) 

25.  WRITE  (6.5200) 

36.  WRITE  (6.1000)  (V(NL.I.l).I  *  1.  ITM1 ) 

37.  C  WRITE  (8.1500)  (V(NL.I.l).I  *  1.  ITM1) 

33.  C  WRITE  (8.1500)  ( UR ( I ) . I  *  1.  ITM1 > 

39.  RETURN 

40.  4  CONTINUE 

41.  C*. ....... .SEARCH  FOR  MAXIMUM  INTENSITY  AT  EACH  GRID  POINT  AND  RECORD 

42.  C. ............ .THE  TIME  AT  WHICH  IT  OCCURS 

43.  IF  (IT.GT.IFRQ)  GO  TO  110 

44.  DO  100  1  =  1.  NT  GP .  2 

45.  EMAX ( I )  =  0. 

46.  100  CONTINUE 

47.  110  CONTINUE 

48.  00  120  1  =  1.  NT  GP.  2 

49.  ABSV  =  ABSi  U( I )  ) 

50.  IF  ( ABSV.LT.EMAXI I ) )  GO  TO  120 

51.  EMAX(I)  =  A8SV 

62.  IT  I  ME  (  I  )  *  IT 

53.  120  CONTINUE 

54.  IF  ( IT.LT. ITM)  RETURN 

55.  WRITE  (6.1500)  (EMAX(I).I  =  1.  NT  G° ,  2) 

56.  WRITE  (6.5000) 

57.  WRITE  (6.9000)  (ITIME(I).I  =  1.  NT  GP .  2) 

33.  C  WRITE  (8.1500)  (EMAX(I).I  =  1.  NTGP,  2) 

59.  C  WRITE  (8.9000)  (ITIME(I).t  =  1.  NTGP.  2) 

eo.  return 


61  .  60  CONTINUE 

62*  IP  (X  —  6)  5  •  o*  6 

63.  5  CONTINUE 

64.  C . ...••PRINT  FIElD  AT  POSITION  X  =  IT  *  H  FOR  TIMES 

65.  C . ....TAJ  s  IT  *  H.  (IT  ♦  IFRC  *  2)  *  H,  ... 

66.  IMAX  =  1  +  ( ITM  -  (IT  /  2)  )  /  IFRO 

67.  WRITE  ( 6 «  30  0 0 )  IT 

68.  WRITE  (6*10001  (U(I).I  s  l.  IMAX) 

69.  C  WRITE  (3*6000)  IT 

70.  C  WRITE  (8*1500)  (U(I).I  *  1*  IMAX) 

71.  RETURN 

72.  6  CONTINUE 

73.  WRITE  (6*4000) 

74.  RETURN 

75.  1000  FORMAT  (1X.9E10.3) 

76.  1500  FORMAT  (8E10.3) 

77.  2C00  FORMAT  <6I5*/*6I5) 

78.  3C00  FORMAT  (6E12.5) 

79.  4000  FORMAT  (///♦ 1 X . 3HF INI SHED ) 

80.  5C00  FORMAT  { / // . IX . 1 4H INC I DENT  FIELD*/) 

81.  5100  FORMAT  ( /// , 1 X • 1 5HREFLECTED  FIELD,/) 

82.  5200  FORMAT  ( /// . 1 X *  I THTR ANSMI TTED  FIELD*/) 

£3.  6000  FORMAT  (  I  10) 

84.  7000  FORMAT  (///* 1 X « 29H INTERNAL  FIELD  AT  TIME  TAU  *  *110. SH  *  2  *  H,/) 

85.  8000  FORMAT  ( /// • IX  * 21 MFI ELD  AT  POSITION  X  **ilO,4N  w  H*/> 

86.  9000  FORMAT  (8110) 

87.  END 
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XII.  PLOTS  FROM  SAMPLE  PROGRAMS 


Internal,  reflected,  and  transmitted  fields  for  mediums  with  nonzero  and  zero 
conductivity.  Incident  field  is  1000-MHz  pulse  of  duration  20  ns,  impinqing 
on  the  medium  from  the  left  at  time  t=0. 


t=6  ns 
t =21  ns 


Nonzero  conductivity  Zero  conductivity 


Figure  A-4 
"  A-5 


Figure  A-10 
"  A-ll 


Comparison  of  transient  and  steady  state  intensities  in  mediums  with  nonzero 
and  zero  conductivity.  Incident  field  is  1000-MHz  pulse. 


Nonzero  conductivity 
Figure  A-6 

Internal  fields,  nonzero  and  zero  conductivity. 

Nonzero  conductivity 


Zero  conductivity 
Figure  A-12 


Zero  conductivity 


x=. 15625 

Figure  A-7 

Figure 

A- 13 

x=.25 

"  A-8 

tl 

A-14 

x=. 71875 

“  A-9 

<1 

A-15 

Reflected  and  internal  fields,  zero  conductivity,  spike  incident  field. 


t= . 181 

ns 

Figure  A- 16 

t=.361 

ns 

"  A-17 

t=.542 

ns 

"  A-18 

t=.723 

ns 

"  A-19 
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I  .  JU  U  •  UU  U  •  SO  1.00  l.SO  2.00 

NORMALIZED  OEPTH 
Figure  B-4.  Internal,  reflected,  and  transmitted  fields  for  a  medium  with 
nonzero  conductivity,  time  t=6  ns.  Incident  field  is  1000-MHz 
pulse  of  20-ns  duration,  impinging  on  the  medium  from  the  left  at 
time  t=0. 
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Figure  B-5.  Internal,  reflected,  and  transmitted  fields  for  a  medium  with 
nonzero  conductivity,  time  t=21  ns.  Incident  field  is  1000-MHz 
pulse  of  20-ns  duration,  Impinqlng  on  the  medium  from  the  left  at 
time  t*0. 
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UO  1.20  2.00  2. 60  3. SO 
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Figure  B-6.  Comparison  of  transient  and  steady  state  intensities  in  a  medium 
with  nonzero  conductivity.  Incident  field  is  1000-MHz  pulse. 
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Figure  B-7.  Internal  field  (x=.!S625),  nonzero  conductivity. 


NORMALIZED  TIME 

Internal  field  (x=. 71875),  nonzero  conduct! 
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NORMALIZED  DEPTH 

Figure  B-10.  Internal,  reflected,  and  transmitted  fields  for  a  medium  with 
zero  conductivity,  time  t=6  ns.  Incident  field  is  1000-MHz 
pulse  of  20-ns  duration,  impinging  on  the  medium  from  the  left 
at  time  t=0. 


-°nWISo  OE^'ffi  100  *■ *  ” 


Figure  B-ll.  Internal,  reflected,  and  transmitted  fields  for  a  medium  with 
zero  conductivity,  time  t*21  ns.  Incident  field  Is  1000-MHz 
pulse  of  20-ns  duration.  Impinging  on  the  medium  from  the  left 
at  time  t*0. 
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Figure  B-12.  Comparison  of  transient  and  steady  state  Intensities  In  a  medium 
with  zero  conductivity.  Incident  field  is  1000-MHz  pulse. 
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Figure  B-15.  Interna)  field  (x=.7IR75),  zero  conductivity. 
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Figure 


Figure 


B-16.  Reflected  and  internal  fields  (t=.181  ns),  zero  conductivity 
spike  incident  field. 
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B-17.  Reflected  and  internal  fields  (t= .361  ns),  zero  conductivity, 
spike  Incident  field. 
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Figure  B-18 .  Reflected  and  internal  fields  (t=.542  ns),  zero  conductivity, 
spike  incident  field. 
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Figure  B-19.  Reflected  and  internal  fields  ( t= .723  ns),  zero  conductivity, 
spike  incident  field. 


